T10N  PAGE 


form  AOfrov+4 
OUt  Ma  0/04-0 IU 


AD-A258  693 


if>«  i»itnw«  m  i«iorw,»no«.  W>«  nytrmrq  cm  »w<  inm  «r  m  ccmr  mm  rf  cm 

i  wwm  D«m#m  N»  *»*«««  Owmmi  * -T  i  ■.  ii  it  mtmw 

M  lo«^n.  i ww^i,  et  leni. 


1.  AGINCT  USI  ONCT  (U«v«  Win*) 


|'-",?}Sv“»,  1992 


J.  MPORT 


ssi  'mmsr 


4.  T1TLI  AND  JUITTPJ 

5.  FUNDING  NUMIIU  f  j 

Theoretical  and  STM  Studies  of  the  Electronic  Structure  of 
Metal/Semiconductor/Hydrogen  Systems 

A&gzpzs  ^ ^ 

?! 'OlUr 

i.  AUTHOKIJJ 

Prof.  Walter  C.  Ermler 

7.  ORGANIZATION  NAMI(S)  AND  ADORIIStlS) 

Ai 

4.  FIRFORJMNG.OROAJUZAnON 

FI  PORT  NUMCIR 

no:;  -r 

y  ns 

Stevens  Institute  of  Technology 

Hoboken,  NJ  0703 1 

•  ->4 

5-27081 

1.  SPONSORING /MO  NIIOIUN  Q  AGINCT  NAMI(S)  ANO  ADOAI5SUS) 

10.  SPONIOWte/MONTTORtNa. 
AUMCT  MPOCT.  NUMUX . 

AJOS’R/*& 

■SUh.  410 

n,  'DC  20332 

^30I/A7 

11.  SUrrUMINTART  NOUS 

•  • 

• 

J 

DTIC 

C*t  P/»>pr"  IMBA. 

ut.  otsTxminoM/AvAiutuurr  jtatimint 


(_;■  f\\  Wt\^  Ve  c\ 


EC  2  91992 


1J.  AiSTRACT  (Minmum JOOwordt) 

The  structure  of  atomic  clusters  is  investigated  using  methods  of  ab  initio  quantum  chemistry  in 
conjunction  with  scanning  tunneling  microscopy  (STM).  Gold  sols  are  imaged  via  STM  and  their  sizes 
characterized.  Atomic  resolution  reveals  facets,  steps,  and  reconstructed  surfaces  of  the  sols.  STM  images  of 
graphite  surfaces  are  characterized  by  two  mechanisms:  a  twisted  top  layer  configuration  and  a  graphite-flake- 
contaminated  tip.  A  crystalline  tip  model  reproduces  anomalous  long  range  periodicity,  attributed  to  defect- 
mediated  tip-subsrate  convolution.  Ab  intio  calculations  of  elecronic  states  of  a  135  atom  Be  cluster  show  bulk 
behavior  for  many  properties  centered  at  the  middle  of  the  cluster.  The  density  of  states  agrees  well  with  band 
theory  models.  Full  spin-orbit  Cl  calculations  of  the  electronic  spectrum  and  potential  energy  curves  of  LiBe, 
including  Rydberg  slates,  are  definitive  and  accurate.  A  procedure  for  the  incorporation  of  core/valence 
polarization  effects  in  all-valence  electron  calculations  is  presented  for  use  with  relativistic  effective  core 
potentials.  The  method  allows  the  treatment  of  large  clusters  at  the  ab  initio  level  including  procedures  for  the 
inclusion  of  electron  correlation. 
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3.  Research  Objectives  and  Status  of  Research 

Grant  AFOSR-90NP123 

Theoretical  and  STM  Studies  of  the  Electronic  Structure  of 
Metal/Semiconductor/Hydrogen  Systems 

Walter  C.  Ermler,  Professor  of  Chemistry  and  Physics 
Department  of  Chemistry  and  Chemical  Engineering 
Stevens  Institute  of  Technology 
Hoboken,  New  Jersey  07030 
(201)  216-5520 


Final  Report 

A  research  project  to  study  the  physical  properties  of  gold  sols  using  the  scanning 
tunneling  microscope  (STM)  was  carried  out  by  a  graduate  student  receiving  support 
through  this  grant  under  the  direction  of  the  PI.  The  experimental  work  was  completed 
within  the  duration  of  the  grant  and  the  student  is  finalizing  his  disseration  and  is  planning  to 
defend  it  this  winter.  A  copy  of  the  dissertation  will  be  forwarded  to  the  AFOSR  upon  its 
final  approval.  The  research  involved  the  deposition  of  colloidal  gold  onto  highly  ordered 
pyrolytic  graphite  and  imaging  using  the  STM.  Two  preparations  were  used  involving  the 
citrate  reduction  of  gold  tetrachloroaurate  yielding  highly  monodispersed  sols.  One  sol 
synthesized  at  373  K  was  found  to  have  a  mean  particle  size  of  16.64  nm  with  a  spread  of 
3.84  nm  from  an  analysis  of  194  particles.  The  mean  diameter-to-height  ratio  was 
detertmined  to  be  8.0  ±  4.6  nm.  A  gold  sol  synthesized  at  346  K  yielded  a  mean  particle 
diameter  of  17.55  ±  4.30  nm  and  diameter-to-height  ratio  of  8.4  ±  2.2.  Examination  of  the 
particle  surfaces  using  the  STM  revealed  atomic  resolution  of  facets,  steps  and  reconstructed 
surfaces.  The  lattice  constants  agreed  well  with  independently  determined  bulk  values. 

Four  additional  related  projects  were  partially  supported  by  this  grant.  The  first 
involved  the  study  of  the  nature  of  STM  images  on  graphite  surfaces  ana  resulted  in  a 
publication  in  The  Journal  of  Physical  Chemistry  on  interpretation  of  apparent  long  range 
order.  Computer  simulations  were  compared  with  actual  STM  data  to  lend  credence  to  the 
postulated  model.  This  paper  is  given  in  the  following  section. 
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The  second  project  dealt  with  the  interpretation  of  large  metal  clusters  in  terms  of  bulk 
properties.  A  cluster  of  135  beryllium  atoms  was  treated  using  effective  core  potentials  (EP) 
to  replace  the  Be  Is  electrons  and  ab  initio  Hartree-Fock-Roothaan  SCF-MO  theory  to 
describe  the  remaining  270  valence  electrons.  Ionization  potentials,  binding  energies,  orbital 
energies,  and  one-electron  properties  were  calculated  for  a  range  of  electronic  states.  The 
ground  state  valence  wave  function  was  also  used  to  calcuate  the  density  of  states  for 
comparison  to  band  structure  calculations.  The  convergence  of  properties  to  bulk  values  was 
determined  based  on  values  of  properties  centered  at  the  atom  in  the  middle  of  the  cluster. 

The  complete  study,  published  in  the  International  Journal  of  Quantum  Chemistry,  is 
discussed  under  Publication  2  below. 


The  use  of  EPs  to  replace  the  Is  core  electrons  of  Li  and  Be  was  the  basis  of  full 
spin-orbit  configuration  interaction  calculations  on  the  prototype  simplest  heteronuclear  metal 
cluster  LiBe  and  three  lowest-lying  cations.  The  complete  electronic  spectrum  including  the 
3s,  3p,  3d  and  4s  Rydberg  states  of  Li  was  calculated  and  compared  to  available  experimental 
results.  Excitation  energies  agree  well  with  experiment  and  the  predicted  spectrum  and 
potential  energy  curves  are  definitive.  The  results  were  published  in  The  Journal  of  Chemical 
Physics  and  are  given  below  under  Publication  3. 

The  LiBe  study  was  an  initial  example  of  a  method  devised  for  treating  few-electron 
systems  at  a  high  level  of  theory.  An  extension  of  the  method  to  include  core/valence 
polarization  is  shown  below  under  Publication  4  and  will  be  published  in  the  Proceedings  of 
a  NATO/ASI  Institute.  A  method  for  deriving  a  core/valence  polarization  operator  that 
includes  spin-orbit  and  other  relativistic  effects  is  presented.  This  operator  can  be  used 
together  with  relativistic  effective  core  potentials  to  reduce  the  number  of  valence  electrons 
that  must  be  explicitly  considered  to  the  smallest  possible  number. 

The  four  projects  were  completed  during  the  period  of  the  one-year  AFOSR  award. 
They  all  contributed  to  a  deeper  understanding  of  the  fundamental  nature  of  interatomic 
interactions.  The  work  was  published  and  the  benefits  to  the  students  involved  in  the 


projects  were  many  fold.  Results  were  also  presented  at  meetings  and  in  departmental . . 
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On  the  Interpretation  of  Scanning  Tunneling  Microscope 
Images  Showing  Anomalous  Periodic  Structures 

Makoto  Sawamura,'  John  F.  Womelsdorf,  Walter  C.  Ermler. 

Department  of  Chemistry  and  Chemical  Engineering, 
Department  of  Physics  and  Engineering  Physics, 
Stevens  Institute  of  Technology, 

Hoboken,  New  Jersey  07030 

Abstract 

Scanning  tunneling  microscope  images  of  graphite  surfaces 
containing  anomalous  long-range  structures  are  shown  to 
originate  from  two  distinct  mechanisms.  A  twisted  top  layer 
graphite  configuration  is  shown  to  produce  super  structures 
with  divergent  periods  from  101  A  to  760  k  and  decaying 
amplitudes.  A  graphite-flake-contaminated  tip  is  shown  to 
yield  long  range  ordered  arrays  with  constant  periodicity, 
constant  amplitudes,  and  abrupt  domain  termination.  These 
arrays  are  consistently  reproduced  using  a  theoretical 
crystalline  tip  model.  Hexagonal  closed-packed  patterns 
with  periods  ranging  from  48  to  220  k  are  experimentally 
observed  and  theoretically  simulated  using  this  model.  The 
origin  of  such  images  is  discussed  in  detail  and  a  class  of 
the  anomalous  long-range  periodicities  observed  are 
attributed  to  a  defect-mediated  tip-substrate  convolution 
phenomenon . 
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Introduction 

Since  its  introduction,^-  the  scanning  tunneling 
microscope  (STM)  has  provided  direct  images  of  surfaces  of 
various  conductive  materials . A  compelling  feature  of 
STM  is  its  ability  to  investigate  local  surface  features  of 
various  substrates  with  atomic  resolution.  However,  the 
detailed  interpretation  of  some  classes  of  images  is  not 
sufficiently  understood.  For  example.,  anomalous  long-range 
periodicity  is  occasionally  reported  to  be  observed  on 
highly  oriented  pyrolitic  graphite  (HOPG)^"®  and  gold. 3 
Although  the  existing  theories  of  the  STM  measurement  do  not 
completely  explain  these  phenomena,  there  are  several 
qualitative  explanations  for  the  long-range  periodicity.6'7 

Mizes  et  al.  first  invoked  the  possibility  of  using 
multiple  tips  for  an  STM  model  to  explain  various  anomalous 
surface  images  of  HOPG.^  Albrecht  et  al.  attributed 
anomalous  long-range  periodicity  on  a  graphite  surface  to  an 
artifact  due  to  isolated  multiple  tips  simultaneously 
scanning  two  grains  in  which  crystal  axes  are  twisted  with 
respect  to  one  another.6  On  the  other  hand,  similar 
anomalous  long-range  periodicity  is  also  believed  to  result 
from  a  twisted  top  layer  configuration  resulting  in  an  STM 
image  described  as  a  moire  fringe  pattern.7  Such  anomalous 
images  are  described  by  constant  periodicity  and  amplitude 
with  abrupt  domain  termination. 


The  origin  of  anomalous  long-range  structures  observed 
on  HOPG  by  STM  is  addressed  in  this  work.  It  is  clearly 
shown  that  a  twisted  top  layer  HOPG  configuration  explains 
divergent  long  range  structures  but  is  inconsistent  with 
ordered  arrays  observed  on  HOPG.  Based  on  a  fundamentally 
stable  model  of  a  polyatomic  crystalline  STM  tip 
configuration,  which  is  logically  extended  to  graphite- 
contaminated  tips,  STM  images  of  such  ordered  arrays  are 
correctly  described.  This  model  shows  that  such  images  are 
an  artifact  of  the  measurement  process  and  not  a  property  of 
the  surface.  The  computational  simulation  based  on  the 
present  model  is  found  to  be  consistent  with  the 
experimentally  observed  anomalous  long-range  order  in  STM 
images . 

Theory 

Many  of  the  existing  theories  of  STM  have  assumed  that 
an  ideal  tip  is  comprised  of  a  single  atom  at  the  position 
nearest  to  the  surface,  although  the  stability  of  the 
single  atomic  tip  has  been  questioned. in  recent  studies 
by  transition  electron  microscopy,  crystalline  structures 
were  observed  on  the  surfaces  of  STM  tips.^'*  The  polyatomic 
STM  tip  model  given  here  derives  from  currently  accepted 
fundamental  assumptions^  containing  the  following 
components.  First,  an  STM  tip  is  assumed  to  possess  a 
finite  crystal  structure,  which  is  inherently  more  stable 
than  a  single  atom  tip  configuration.  Second,  the  net 
tunneling  current  is  defined  as  the  summation  of  currents 


due  to  tunneling  from  each  atom  on  the  tip  to  every  atom  of 
the  sample  substrate.  The  intensity  of  the  current  between 
an  atom  on  the  sample  and  an  atom  on  the  tip  is  described  by 
the  Tersoff  and  Hamann  relation  in  the  low  bias  limit.1® 

The  net  tunneling  current  from  an  STM  tip  positioned  at 
the  coordinate  ( l,k )  to  the  surface  can  be  defined  by: 


I ijU'k)  -  A  exp[-a{  (xijj-xIjc)2+{yij-y2^)2+(2ij-zj*)2}1/2] , 


where  I ±j(l,k)  is  the  tunneling  current  between  an  atom 
located  at  on  the  surface  and  an  atom  at 

(xj^,  yik1  zlk)  on  the  STM  tip,  and  a  depends  on  the 

characteristics  of  the  materials  of  both  the  tip  and 
surface.  The  net  current  is  a  summation  of  all  tunneling 
occurrences  resulting  in  the  superimposed  image.  In  the 
case  where  the  sample  substrate  has  a  low  workfunction 

region,  for  example  a  step  edge,  17  it  has  been  shown  that 
the  defect  may  act  as  a  pseudo-adatom  and  be  scanned  by  the 
STM  tip  surface.1®  In  this  instance  the  atoms  on  the  tip 
will  periodically  scan  the  point  defect  producing  a  tip 

self-image.  The  analysis  presented  in  this  paper  will 

assume  that  graphite  step  edges  constitute  low  barrier 
height  regions.  Consequently,  it  is  therefore  not 

surprising  that  in  many  cases  long-range  periodicity  has 
been  observed  adjacent  to  a  surface  defect.^-® 


Results 


Figure  1  displays  an  image  of  a  long-range  super 
structure  observed  on  HOPG  in  air  which  is  attributed  to  a 
twisted  top  layer  HOPG  configuration.  The  bias  voltage  was 
50  mV  and  the  tunneling  current  1.0  nA.  The  data  was 
obtained  using  a  PtQ.8Ir0.2  tip.  The  observed  periodicities 
range  from  101  A  to  760  A  smoothly  diverging  along  the  sides 
of  the  triangle  domain.  This  image  clearly  describes  the 
divergent  nature  of  the  anomalous  image  which  is  consistent 
with  a  twisted  layer  mechanism  based  on  the  discussion 
below. 

Figure  2a  displays  an  STM  image  of  a  long-range 
periodic  pattern  on  a  freshly  cleaved  HOPG  surface  probed  in 
air  with  a  Ptg  giro  2  tip*  The  bias  voltage  was  150  mV  and 
the  tunneling  current  set  at  1.0  nA.  The  STM  was  operated 
in  constant  height  mode.  A  hexagonal  close-packed  pattern 
with  long-range  periodicity  of  60  +/-  1  A  is  apparent.  It 
is  noted  that  atomic  resolution  can  be  observed  throughout 
the  image.  The  large  hexagons  and  atomic  hexagons  are 
rotated  with  respect  to  one  another  by  approximately  30 
degrees.  This  image  displays  a  portion  of  the  surface  which 
contains  a  long-range  periodic  domain  lying  along  a  step 
edge. 


A  3-dimensional  simulation  for  the  experimental  image 
of  Figure  2a  based  on  a  tip-substrate  convolution  is  shown 
in  Figure  2b.  The  sample  substrate  was  taken  as  a  hexagonal 
close-packed  surface  with  an  atomic  distance  of  2.45  A 


corresponding  to  the  HO  PG  (0001)  surface  in  keeping  with 
known  a  and  (5  inequivalences  on  graphite.19  Both  the  sample 
substrate  and  tip  are  comprised  of  25  x  35  =  875  atoms.  In 
this  simulation,  a  graphite  (0001)  surface  was  assumed  for 
both  the  tip  and  the  sample  substrate  and  includes  P  atoms 
only.  The  substrate  contains  a  low  workfunction  defect  in 
the  scanned  region.  The  coordinate  axes  of  the  tip  and 
sample  are  rotated  relative  to  one  another  by  2.3  degrees. 
A  large  periodicity  is  clearly  observed  in  Figure  2b  having 
a  period  of  60  A.  This  corresponds  to  the  experimentally 
obtained  result  of  60  +/-  1  A  shown  in  Figures  2a. 

It  is  clearly  plausible  to  envisage  a  graphite  surface 
for  an  STM  tip  instead  of  a  metal  surface  since  the  tips  are 
known  to  pick  up  graphite  flakes.2*^21  After  scanning  a 
graphite  surface  for  a  sufficient  period  of  time  the 
graphite  is  observed  to  grow  hair-like  on  the  STM  tip. 
Since  these  graphite  hairs  must  be  in  contact  with  the 
surface  substrate  as  they  grow,  the  graphite  surface 
contributing  to  the  net  tunneling  current  may  be  nearly 
parallel  to  the  sample  substrate.  This  is  conceptually 
similar  to  a  paintbrush  that  is  placed  in  close  contact  with 
a  planar  surface  which  will  relax  the  bristles  to  a  nearly 
coplanar  configuration  in  order  to  reduce  the  repulsive 
forces.  Such  a  graphite  flake  results  in  a  large  periodic 
tip  surface  which  explains  the  large  ordered  super 
structure.  The  super  periodicity  s  for  the  graphite 
substrate  and  tip  is  described  by  3  -  a/2sin(0/2)  where  a  is 
the  graphite  distance  and  9  the  rotational  angle  of  the 


crystal  axis  of  the  tip  with  respect  to  the  substrate 
crystal  axis.  A  similar  expression  for  the  graphite  super 
periodicity  can  be  found  in  the  work  of  Kuwabara  et  al."^ 
However,  the  mechanism  presented  here  is  completely  distinct 
from  their  model.  The  similarity  between  both  formulations 
lies  only  in  their  geometric  descriptions. 

Discussion 

A  twisted  top  layer  model  of  HOPG  may  be  divided  into 
two  cases .  The  first  case  is  that  of  a  large  flake  which  is 
rotated  at  some  angle  relative  to  the  basal  plane  of  the 
second  layer.  The  second  case  corresponds  to  a  region  of 
the  top  layer  basal  plane  which  has  been  rotated  in  a 
similar  fashion  but  represents  a  small  domain  of  an 
otherwise  commensurate  layer.  The  former  configuration  is 
dismissed  since  it  requires  that  the  domain  of  the  anomalous 
structure  be  contained  within  the  step  edges  of  the  flake. 
Such  structures  have  not  been  reported.  The  latter, 
however,  does  not  demand  a  step  edge  perimeter  although  the 
presence  of  such  a  step  edge  may  indeed  be  present. 

Although  it  has  been  reported  that  a  twisted  top  layer 
configuration  results  in  long  range  periodic  STM  images, 
such  reports  show  distinct  boundaries  between  the  super 
structure  and  the  normal  graphite  lattice.^  There  is 
usually  a  defect  associated  with  one  boundary  with  the 
others  terminating  abruptly.  This  suggests  that  the 
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rotation,  which  has  been  proposed  as  an  explanation  for  such 
features,  also  terminates  abruptly.  This  logically  results 
in  a  tearing  of  the  top  layer,  which  is  not  experimentally 
observed. 

However,  if  a  top  layer  basal  plane  contains  two 
domains,  one  twisted  with  respect  to  the  underlying  layer 
and  the  other  left  unperturbed,  then  the  maximum  rotation 
present  in  one  domain  must  be  damped  to  zero  rotation  in  the 
other  domain.  This  transition  will  avoid  cleaving  the  basal 
plane  which  contains  the  two  domains,  but  requires  that  the 
angle  of  rotation  changes  over  the  transition  or  damped 
region.  Therefore,  the  predicted  periodicity  cannot  be 
constant  over  this  region.  Furthermore,  the  effective  angle 
of  rotation  over  this  transition  region  is  gradually  reduced 
as  the  unperturbed  domain  is  approached.  As  the  large 
periodicity  increases,  the  rotational  angle  decreases 
according  to  S  =*  a/2sin(0/2)  as  proposed  by  Kuwabara  et  al.^ 
Since  the  wavelength  of  the  long-range  periodicity 
approaches  infinity  at  infinitesimal  rotation  angles,  the 
boundary  shared  by  the  two  domains  is  described  by  the  limit 
of  the  diverging  superstructure.  Figure  1  displays  just 
this  behavior.  The  superstructure  is  described  by  an  array 
whose  long-range  period  diverges  as  the  normal  graphite 
lattice  is  recovered  accompanied  by  the  amplitude  of  the 
array  damping  to  zero.  Therefore,  Figure  1  is  consistent 
with  a  twisted  top  layer  model. 


It  was  shown  in  our  previous  study  that  if  the  tip  and 
sample  have  different  nearly  parallel  crystalline 
structures,  and  there  is  no  low  workfunction  defect  on  the 
surface  corresponding  to  the  smaller  atomic  lattice,  the 
observed  image  is  a  periodic  pattern  of  the  surface 
containing  the  smaller  periodicity.1®  This  result  may  be 
apparent  if  translationally  symmetric  properties  of  the 
system  are  considered  in  the  scanning  process.  Albrecht  et 
al.  employed  isolated  multiple  tips  scanning  two  grains 
twisted  relative  to  one  another.®  However,  images  having 
two  domains  exhibiting  both  normal  atomic  resolution  and 
long-range  periodicities  along  the  grain  boundaries  cannot 
be  explained  using  their  model  {Figure  2  a) .  First  of  all, 
it  is  unlikely  that  two  minitips  resolving  atomic 

periodicity  contribute  to  the  net  tunneling  current  on  a 
nearly  equal  basis.  Second,  it  is  unclear  how  distinct  well 
defined  boundaries  will  emerge  if  the  grain  boundary  is 
irregular  as  was  seen  in  Ref  6.  Finally,  this  model  does 
not  address  why  such  long-range  structures  are  not  observed 
on  other  surfaces  which  exhibit  grain  boundaries. 

Long-range  order  over  three  different  steps,  where 
these  steps  produced  identical  long-range  spacings,  is  shown 
in  Figure  3.  The  surface  was  scanned  in  air  with  a 
PtQgIro.2  tip  under  20.1  mV  bias  voltage  and  1.0  nA 
tunneling  current  in  constant  height  mode.  The  period  is  52 
+  /“  2  A.  This  image  is  attributed  to  a  tip-substrate 

convolution  resulting  from  a  polyatomic  STM  tip  scanning  an 
HOPG  substrate  which  contains  a  low  workfunction  defect . 


Kuwabara  et  al.,  however,  attributed  the  observed  long-range 
periodicity  to  the  presence  of  a  twisted  top  layer  of 
graphite. ^  if  a  twisted  first  layer  model  is  used  for  the 
interpretation  of  this  image,  it  is  necessary  to  assume  that 
every  step  containing  the  identical  long-range  periodicity 
is  rotated  by  the  same  angle  with  respect  to  the  next  lower 
layer.  Furthermore,  as  discussed  above,  it  is  unclear  how 
distinct  boundaries  and  constant  periodicity  can  be 
explained  by  this  method.  Consequently,  this  phenomenon 
cannot  be  justified  in  general  using  a  twisted  first  layer 
model.  In  addition  independent  surface  techniques  would  not 
detect  such  long-range  order  since  it  results  as  an  artifact 
of  the  STM  measurement  process. 

Conclusions 


Anomalous  long-range  periodic  structures  observed  on 
HOPG  may  be  divided  into  two  classes.  The  first  class  of 
images  results  from  a  twisted  top  layer  configuration.  The 
nature  of  these  images  is  in  stark  contrast  with  those 

images*:  derived'  from  twisted'  top 
layer  HOPG  structures  are  identified  by  a  long-range  super 
structure  whose  periodicity  diverges  with  a  decaying 
amplitude  as  the  non  twisted  domain  is  approached. 


The  second  class  of  images  result  from  a  crystalline  or 
graphite  contaminated  tip  that  encounters  a  low  workfunction 
surface  defect.  These  images  are  fundamentally  different  in 
that  they  have  constant  long-range  periods. 


constant 
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amplitudes,  and  abrupt  domain  boundary  termination.  Such 
images  are  indicative  of  extensive  tip  contamination.  It 
should  be  noted  that  all  the  images  of  this  type  which  have 
been  encountered  in  our  laboratory  may  be  successfully 
interpreted  based  on  this  model. 
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Figure  Captions: 

Figure  2.  A  2500  A  x  2500  A  experimental  image  of  an  HOPG 
substrate  with  bias  voltage  of  50  mV.  The  long-range  period 
diverges  along  the  edge  of  the  triangle  domain,  ranging  from 
101  A  to  760  A. 

Figure  2a.  A  160  A  x  120  A  experimental  image  of  an  HOPG 
substrate,  scanned  in  air  with  bias  voltage  of  150.1  mv. 
The  dark  areas  comprise  a  large  hexagonal  closed-packed 
pattern  with  a  period  of  60  +/-  1  A.  The  atomic  resolution 
is  also  observed. 

Figure  2b.  A  simulation  of  a  graphite  (0001)  surface 
scanned  by  a  graphite  (0001)  surface  comprising  a  tip.  The 
size  of  the  image  is  58.8  A  x  83.3  A.  The  long-range  period 
of  60  A  is  clearly  reproduced.  The  image  is  horizontally 
compressed  by  37.4  %. 

Figure  3.  A  4200  A  x  4200  A  experimental  image  of  an  HOPG 
substrate,  scanned  in  air  with  bias  voltage  of  20.1  mV.  The 
long-range  period  of  52  +/-  2  A  is  observed  in  two  domains 
over  three  steps . 
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Abstract 

Hartrec-Fock-Roolhaan  studies  arc  reported  for  low-lying  electronic  states  of  metallic  beryl¬ 
lium  as  modeled  by  a  moiety  of  135  beryllium  atoms.  The  system  corresponds  to  16  coordination 
shells  of  a  central  Be  with  internuclear  separations  derived  from  the  lattice  constants  of  the  bulk 
metal.  The  calculations  become  tractable  by  use  of  the  full  Du,  symmetry  of  the  system  at  both 
the  integrals  and  self-consistent-field  stages  and  by  employing  ab  initio  effective  potentials  for 
the  Is  electrons  of  each  beryllium  atom.  Ionization  potentials,  binding  energies,  orbital  energies, 
electric  field  gradients,  nuclear-electrostatic  potentials,  diamagnetic  shielding  constants,  second 
moments,  and  Mulliken  populations  are  calculated  for  selected  electronic  states.  The  calculated 
ionization  potential  for  the  lowest  state  agrees  to  within  10%  of  the  experimental  bulk  work  func¬ 
tion.  A  density-of-states  analysis  for  that  state  is  reported  and  compared  with  band  structure 
calculations. 

Introduction 

Bulk  systems  can  be  modeled  employing  first-principle  quantum  chemical 
studies  on  large  fragments  of  the  bulk  provided  a  sufficient  number  of  atoms  are 
included.  Since  ab  initio  quantum  chemical  calculations  on  systems  of  50- 
150  atoms  arc  becoming  increasingly  manageable  with  current  supcrcomputing 
technology  and  systems  of  larger  dimension  arc  likely  to  become  feasible  with 
exploitation  of  massively  parallel  systems,  it  is  of  interest  to  test  this  approach  on 
a  bulk  system. 

In  choosing  a  model  system  for  exploratory  theoretical  study,  the  availability 
of  experimental  results  is  an  important  criterion.  In  addition  to  the  fact  that 
many  of  its  bulk  properties  as  well  as  its  electric  field  gradient  arc  known,  beryl¬ 
lium  has  attracted  a  considerable  amount  of  theoretical  interest.  Previous  large 
fragment,  cluster,  or  band  theory  studies  include  ab  initio  Hartrec-Fock- 
Roothaan  (hfr)  studies  on  fragments  of  the  bulk  metal  of  up  to  87  atoms  [1-5], 
hfr  studies  on  a  cluster  of  55  atoms  [6],  a  Hartree-Fock  band  theory  study  [7], 
local  density  approximation  band  theory  studies  [8,9],  and  orthogonalized  plane 
wave  band  theory  studies  [10.11]. 
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In  this  work  Hartrcc-Fock-Roothaan  calculations  arc  reported  for  selected 
states  of  metallic  beryllium  as  modeled  by  a  moiety  of  135  beryllium  atoms. 
A  multiconfiguration  self-consistent  field  (mcscf)  calculation  is  also  reported  to 
investigate  the  mixing  of  two  low-iying  states.  The  system,  shown  in  Figure  1,  is 
a  fragment  of  bulk  beryllium  metal  and  corresponds  to  16  coordination  shells 
of  atoms  around  a  central  beryllium  atom  (13  shells  if  bulk  beryllium  were  per¬ 
fect  hep).  This  work  is  an  extension  of  previous  studies  on  metal  fragments  of 
13  through  87  beryllium  atoms  [1-5].  Results  are  compared  to  experiment,  the 
earlier  fragment  studies,  and  studies  employing  theoretical  methods  based  on 
band  theory. 


Calculations 

The  iifr  calculations  become  tractable  through  use  of  the  full  Dih  point  group 
symmetry  of  the  system  and  ab  initio  effective  potentials  to  replace  the  Is  core 
electrons  of  each  beryllium  atom.  A  contracted  (3s2p)/[2slp]  Gaussian-type 
function  (gtf)  basis  set  optimized  in  a  previous  study  [1]  has  been  employed  to 
represent  the  valence  electrons.  In  the  integrals  step  of  the  study,  2  x  10“  integrals 
were  computed  requiring  43  h  of  Cray  Y-MP8  supercomputer  cpu  time;  this 
would  be  a  considerable  amount  on  a  one-processor  system,  however,  optimal 
use  of  multiprocessor  systems  such  as  the  Cray  Y-MP8  can  make  such  calcula¬ 
tions  routinely  manageable.  A  comparable  all-electron  study  with  a  [9j5/>]  primi¬ 
tive  gtf  basis  set  and  without  the  use  of  symmetry  or  effective  potentials  would 
require  the  computation  of  2  X  10u  integrals.  Details  of  the  specific  reductions 


Figure  1.  Beu>:  16  coordination  shells  of  a  central  Be  atom. 


in  numbers  of  integrals  with  use  of  symmetry  and  effective  potentials  have  been 
given  in  an  earlier  study  [4]. 

The  entire  set  of  two-electron  integrals  required  4.8  gigabytes  of  disk  space.  In 
subsequent  hfr  and  mcscf  studies,  however,  only  the  subset  of  integrals  required 
for  a  particular  spin  state  was  kept  on-line  for  states  of  that  class.  This  resulted 
in  a  reduction  of  on-line  disk  space  requirements  to  a  minimum  of  3.3  gigabytes. 

The  procedure  followed  initially  for  the  hfr  studies  was  to  guess  a  closed-shell 
electron  configuration  based  on  low-lying  states  found  for  smaller  fragments. 
Molecular  orbital  occupations  were  then  adjusted  in  test  calculations  until  the  or¬ 
bital  energies  of  the  highest  occupied  molecular  orbitals  (homo)  for  all  of  the  6 
irreducible  representations  were  lower  than  that  for  the  lowest  unoccupied  or¬ 
bital  (lumo)  of  the  fragment.  This  state  was  then  converged  (Table  I,  state  9)  to 
10“*  a.u.  for  one-  and  two-electron  energies.  The  resulting  wave  function  was 
then  used  as  an  initial  guess  for  an  exhaustive  search  for  lower-lying  states  using 
homo/lumo  analyses  as  a  partial  guide. 

Molecular  orbital  (mo)  occupations  and  calculated  binding  energies,  ionization 
potentials,  and  excitation  energies  for  the  ground  state  and  7  additional  lower- 
lying  states  converged  to  10"’  a.u.  are  shown  in  Table  I.  These  states  correspond 
to  either  closed-shell  or  average-energy-of-configuration  states  in  the  case  of 
open-shell  mo  configurations.  The  average  energy  of  a  single  mo  configuration  is 
defined  as  the  weighted  mean  of  the  energies  of  all  multiplets  for  the  configura¬ 
tion  [12].  In  earlier  studies  on  fragments  of  greater  than  50  atoms,  it  was  ob¬ 
served  that  calculated  properties,  populations,  and  other  quantities  are  largely 
insensitive  to  spin  coupling  for  a  given  electron  configuration  [2]. 

The  hfr  studies  typically  required  between  300  and  400  iterations  to  converge 
each  state  such  that  properties  were  constant  to  the  significant  figures  reported 
(10*’  a.u.  in  this  case).  Each  iteration  required  136  seconds  of  Cray  Y-MP  cpu 
time.  The  hfr  and  mcscf  studies  for  the  various  states  and  ions  were  all  carried 
out  using  the  "once-computed"  integrals  data  set  described  above. 

While  alternative  methods  such  as  direct  scf  [13]  can  be  performed  with 
smaller  disk  space  requirements,  they  were  not  considered  in  this  study  because 
it  was  planned  to  investigate  a  large  number  of  neutral  and  ion  states.  Recompu¬ 
tation  of  the  integrals  in  every  case  would  be  required  if  such  direct  methods 
were  employed.  While  the  current  approach  requires  the  use  of  more  disk  space 
than  direct  methods,  it  results  in  greater  flexibility  with  regard  to  the  investiga¬ 
tion  of  excited  states  and  ions  with  an  appreciable  overall  cpu  time  savings. 

An  mcscf  study  was  carried  out  to  investigate  the  mixing  of  low-lying  elec¬ 
tronic  states.  In  this  case  the  mo  configurations  of  the  ground  state  and  the  state 
obtained  by  exciting  two  electrons  from  highest  occupied  ai  molecular  orbital  of 
the  ground  state  into  the  lowest  unoccupied  a  J  molecular  orbital  were  employed 
as  configurations  for  the  study  (states  1  and  8  of  Table  I).  Total  energies  and  con¬ 
figuration  interaction  coefficients  are  reported  for  this  study  in  Table  II. 

Properties  analyses  were  also  carried  out  for  low-lying  states  as  reported  in 
Table  I.  Those  properties  included  the  electric  field  gradient,  nuclear-electrostatic 
potential,  and  diamagnetic  shielding  constant  at  the  central  atom,  and  the 
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Table  II.  Multiconfiguration  self-consistent  field  study  on  Be,jj 


mcsct  Study 

Configuration 

interaction 

coefficients 

Total  energy  (a.u.) 

State  I 

State  8 

2-Configuration 

0.9997 

-0.0230 

-135.180602 

Molecular  orbital 

scf  Studies 

configuration* 

Total  energy  (a.u.) 

Electronic  state 

19n| 

Closed-shell  state  1 

0 

2 

—  135.180^26 

Triplet 

1 

1 

-135.150690 

2-e"  average  of  configuration 

1 

1 

-135.150496 

Closed-shell  stale  8 

2 

0 

-135.116544 

‘Closed-shell  mos  (lal)j •  •  •  (Ini)1  •  •  •  (701),  (la;)1  •  •  ■  (15 a])\  (le')4 •  •  •  (25e*)4.  (10'~ •  (19e-)‘. 

CI«r)*"*(6eT)s. 


quadruple  and  second  moments  of  charge  for  the  entire  fragment.  Of  these,  the 
calculated  electric  field  gradients  and  quadrupole  moments  for  the  low-lying 
states  reported  in  Table  I  are  shown  in  Figure  _. 
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Figure  2.  Calculated  electric  field  gradients  and  quadrupole  moments  for  low- 
lying  states  of  Bern  within  1.75  cV  of  thc\ground  state  as  reported  in  Table  I.  Prop¬ 
erties  calculated  for  states  5  and  6  are  superimposed. 
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Mulliken  population  analyses  [14]  were  also  carried  out  for  the  low-lying  states 
reported  in  Table  I.  Net  charges  and  gross  atomic  s  and  p  orbital  populations 
were  calculated  for  each  symmetry-distinct  group  of  atoms.  The  average  net 
charges  from  the  analyses  for  the  symmetry-distinct  groups  of  atoms  are  shown 
in  Figure  3. 

Valence  electronic  density  distribution  maps  have  also  been  calculated  for 
the  ground  state.  Contour  levels  were  calculated  starting  at  6.25  x  1CT5 
electrons/a.u.3  and  doubled  thereafter  to  a  maximum  of  65.536  electrons/a.u.3 
The  density  distribution  maps  are  shown  in  Figure  4. 

A  density  of  states  (dos)  analysis  was  also  carried  out  for  the  ground  state  and 
is  shown  in  Figure  5  [15].  To  obtain  the  density  of  states,  a  point  corresponding 
to  the  number  of  states  in  a  given  orb:'al  energy  increment  was  plotted  in  the 
center  of  the  associated  increment.  The  dos  was  then  drawn  directly,  connecting 
the  plotted  points  as  shown  in  Figure  5.  It  was  decided  not  to  fit  the  curves  in 
order  to  avoid  introducing  artifacts. 

The  one-electron  energies  for  the  occupied  orbitals  in  the  dos  were  obtained 
from  the  ground  state  wave  function  while  two  sets  of  energies  have  been  plotted 
for  the  virtual  space.  One  set  of  energies  corresponds  to  that  obtained  from  the 
canonical  Hartree-Fock  (hf)  orbitals  while  a  second  set  has  been  obtained  from 
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Figure  3.  Average  net  atomic  charges  for  the  20  symmetry  equivalent  groups  of 
atoms  as  a  function  of  the  radii  of  the  coordination  shells  [4]  containing  the  respec¬ 
tive  groups.  Sharp  corners  arc  artifacts  of  the  relation  of  groups  with  the  shell 
structure.  Vertical  segments  are  due  to  groups  of  atoms  that  are  s;  mmclry  equiva¬ 
lent  in  the  bulk  but  noncquivalcnl  in  the  fragment. 
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Figure  4.  Contours  of  the  molecular  valence  charge  density  in  the  coordinate 
planes  for  Bc,}s.  Contour  levels  begin  at  6.25  x  10’’  e/(a.u.)s  and  arc  doubled  for 
each  new  contour,  (a)  X Y  plane,  (b)  XZ  plane,  (c)  YZ  plane. 


an  improved  virtual  orbital  calculation  [16,17].  This  was  required  because  the 
Hartrec-Fock  unoccupied  mos  represent  the  system  under  study  plus  one  addi¬ 
tional  electron,  i.e.,  the  negative  ion  [18].  Improved  virtual  orbitals  (ivo)  that 
have  been  corrected  so  that  they  are  consistent  with  the  corresponding  neutral 
system  were  obtained  by  performing  a  single  iteration  self-consistent  field  (scf) 
calculation  on  the  ionized  fragment  with  the  neutral  ground  state  scf  wave  func¬ 
tion  as  input.  In  order  to  obtain  an  overall  average  effect,  an  appropriate  fraction 
of  the  ionized  electron  was  removed  from  each  of  the  highest-occupied  molecu¬ 
lar  orbitals  of  the  irreducible  representations  of  the  D3*  point  group. 

In  order  to  determine  the  optimal  energy  increment,  density  of  states  analyses 
resulting  from  use  of  energy  increments  of  0.01,  0.035,  0.050,  0.065,  0.080,  0.095, 
0.11,  0.135,  0.150,  and  0.175  a.u.  were  examined  (at  the  extreme  limits  in  a  dos 
analysis,  the  smallest  energy  increment  has  a  small  number  or  zero  states  while 
the  largest  contains  the  full  set  of  states).  It  was  found  that  an  energy  increment 
of  0.035  a.u.  contains  the  most  dos  data  points  while  still  maintaining  the  gen¬ 
eral  features  of  the  curves. 
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Discussion 

As  can  be  seen  in  Table  I,  of  the  9  states  converged  to  1Q~ 9  a.u.  in  one-  and 
two-electron  energy  in  the  hfr  studies,  3  are  closed-shell  states,  5  are  open-shell 
average-of-configuration  states,  and  1  is  a  triplet  state.  In  addition  to  these 
9  states,  31  additional  states,  converged  to  0.1  a.u.,  were  found  within  4  cV  of  the 
ground  state.  The  additional  states  were  not  converged  further  because  examina¬ 
tion  of  properties  and  populations  analyses  indicated  little  change  from  those  of 
the  9  states  already  studied  in  detail  (Table  I). 

The  results  of  the  present  study  are  consistent  with  results  of  earlier  fragment 
studies  of  greater  than  50  beryllium  atoms  in  that  calculated  properties,  popula¬ 
tions,  and  other  quantities  are  largely  insensitive  to  spin  coupling  for  a  given 
electron  configuration  [2],  For  example,  as  can  be  seen  in  Table  I,  the  calculated 
ionization  potentional  for  the  triplet  state  is  3.48  eV  and  the  corresponding 
average-of-configuration  ionization  potential  is  3.47  cV.  Population  and  proper¬ 
ties  analyses  yield  similar  comparisons.  For  example,  the  electric  field  gradient  at 
the  central  beryllium  atom  is  calculated  to  be  0.0091  a.u.  for  both  states  and  the 
quadrupole  moments  for  the  entire  fragment  are  calculated  to  be  18.18  and 
18.21  a.u.  for  the  triplet  and  average-of-configuration  states,  respectively.  In  re¬ 
gard  to  population  analyses,  the  calculated  net  charges  with  respect  to  the 
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Figure  4.  (Continued) 

symmetry-distinct  groups  of  atoms  have  been  found  to  be  identical  for  the  triplet 
and  average-of-configuration  states. 

Calculated  ionization  potentials  at  the  scf  level  (Ascf)  for  the  four  lowest 
states,  shown  in  Table  I,  are  all  larger  than,  but  within  10%  of,  the  experimental 
bulk  work  function  [19].  The  Koopmans’theorem  ionization  potential  for  the  low¬ 
est  state  is  0.24  eV  larger  than  the  corresponding  Ascf  ionization  potential  and  is 
15%  larger  than  the  experimental  value. 

The  calculated  binding  energy  per  atom  for  the  ground  state  is  31.7  kcal/mol. 
In  comparing  to  calculated  binding  energies  per  atom  for  the  lowest  states  of 
bulk  fragments  of  51,  57,  69,  81,  and  87  atoms,  the  onset  of  convergence  to  the 
bulk  hfr  limit  employing  this  basis  set  is  observed.  Calculated  binding  energies 
for  the  lowest  lying  states  of  bulk  fragments  of  51,  57,  69,  81,  and  87  atoms  are 
25.8  [2],  25.4  [2],  28.2  [4],  28.3  [5],  and  28.5  [5]  kcal/mol,  respectively.  These  cal¬ 
culated  values  yield  an  incremental  increase  of  0.077  kcal/mol  per  atom  as  the 
number  of  atoms  is  increased  from  51  to  87.  This  can  be  compared  to  a  smaller 
incremental  increase  of  0.066  kcal/mol  per  atom  as  the  number  of  atoms  is  in¬ 
creased  from  87  to  135.  If  the  size  of  the  incremental  increase  is  taken  to  de¬ 
crease  by  at  least  0.01  kcal/mol  per  atom  for  every  42  atom  increase  in  fragment 
size,  an  estimate  of  the  maximum  contribution  to  binding  energy  from  further 
fragment  size  increase  is  about  7  kcal/mol  per  atom. 
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Figure  5.  Density  of  states  (dos)  analyses  for  the  ground  state  of  Bca,  for  an  en¬ 
ergy  increment  of  0.035  a.u. 

An  estimate  of  basis  set  effects  on  binding  energy  can  be  obtained  from  earlier 
studies  on  smaller  beryllium  fragments  and  clusters.  The  effect  of  uncontracting 
the  current  p  basis  set  can  be  obtained  by  comparing  earlier  fragment  studies  on 
Beji  and  Be37  [2]  to  studies  by  Pettersson  and  Bauschlicher  on  Bejs  [6]  and  is  es¬ 
timated  to  increase  the  binding  energy  by  about  3%.  Similarly,  the  effect  ob¬ 
tained  by  adding  an  additional  p  primitive  basis  function  is  estimated  to  increase 
the  binding  energy  by  14%  also  based  on  the  same  studies  [6].  Adding  yet  an¬ 
other  p  primitive  basis  function  is  estimated  to  increase  the  binding  energy  by 
4.5%  based  on  comparing  studies  on  Beu  and  Be3j  again  by  Pettersson  and 
Bauschlicher  [6],  The  effect  of  adding  a  d  function  can  be  estimated  to  increase 
the  binding  energy  by  20%  based  on  an  average  of  effects  observed  in  studies  on 
Be4  by  Marino  and  Ermlcr  [20]  and  Bauschlicher,  Bagus,  and  Cox  [21].  Summing 
these  effects  for  the  ground  state  of  Be  us,  a  maximum  basis  set  error  of 
13  keal/mot  per  atom  is  obtained.  It  is  likely,  however,  that  this  estimate  is  too 
large  based  an  observed  decrease  of  55%  in  the  basis  set  effect  obtained  by 
adding  a  primitive  p  function  to  the  basis  set  in  comparing  clusters  of  13  and 
56  atoms  [6], 

We  attribute  the  remaining  24  to  30  keal/mol  per  atom  required  to  attain  the 
experimental  binding  energy  [22]  to  differential  electron  correlation  energy.  This 
is  supported  by  correlation  effects  seen  in  studies  on  Be«  of  as  large  as  76  and  54% 
by  Marino  and  Ermlcr  [20]  and  Bauschlicher,  Bagus,  and  Cox  [21],  respectively. 
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To  investigate  the  mixing  between  low-lying  states,  an  mcscf  study  was  carried 
out  for  the  two  lowest  closed-shell  states.  As  shown  in  Table  II,  these  states  cor¬ 
respond  to  the  ground  state  and  that  obtained  by  a  double  excitation  from  the 
8«2  orbital  of  this  state  into  the  unoccupied  18a  |  orbital.  The  total  energy  for  the 
mcscf  wave  function,  configuration  interaction  coefficients,  hfr  total  energies 
for  the  two  closed-shell  states,  and  the  triplet  and  average-of-configuration  states 
obtained  from  a  single  excitation  involving  the  same  orbitals  are  shown  in 
Table  II.  As  can  be  seen  from  the  configuration  interaction  (ct)  coefficients, 
there  is  very  little  mixing  of  the  two  states.  This  is  no  doubt  due  to  the  small 
size  of  the  exchange  integral  between  these  two  orbitals  both  of  which  are  spread 
out  over  the  entire  cluster.  Accordingly,  the  energy  lowering  in  the  mcscf  study 
is  only  0.0048  eV  relative  to  the  ground  state,  indicating  that  a  more  extensive 
correlation  treatment  is  necessary  to  recover  significant  correlation  effects, 
as  expected. 

The  calculated  average  values  of  the  nuclear-electrostatic  potential  [23]  and 
diamagnetic  shielding  constant  [23]  are  -1.06  and  -33.10  a.u.,  respectively.  For 
both  properties,  only  one  state  deviated  more  than  0.01  a.u.  from  the  average. 
Consistent  with  properties  computed  for  fragments  of  bulk  beryllium  larger  than 
51  atoms  [2-5],  the  computed  values  are  relatively  constant  with  electronic  state. 
This  behavior  is  also  observed  in  the  calculated  values  of  the  square  roots  of  sec¬ 
ond  moments  of  charge  [23],  (rJ)1/2  and  (zJ) 1/J,  which  are  found  to  have  average 
values  of  93.28  and  88.01  a.u.,  respectively,  and  maximum  state-to-state  varia¬ 
tions  of  0.19  and  0.24  a.u.,  respectively.  As  discussed  in  the  earlier  studies,  this 
can  be  interpreted  as  a  sign  of  bulk  behavior  in  that  smaller  fragments  exhibit 
larger  changes  in  properties  with  electronic  state. 

The  calculated  values  of  the  electric  field  gradient  [23]  at  the  central  atom  and 
the  quadrupole  moment  [23]  for  the  entire  fragment  for  the  lowest  eight  con¬ 
verged  states  are  shown  in  Figure  2  as  a  function  of  excitation  energy.  As  dis¬ 
cussed  above,  the  calculated  values  for  the  electric  field  gradient  and  the 
quadrupole  moment  for  states  5  and  6  (triplet  and  average  of  configuration  for 
the  same  molecular  orbital  occupation)  are  identical  or  close  to  identical  for  the 
number  of  significant  figures  reported.  This  results  in  the  points  for  these  states 
in  the  figure  being  superimposed.  Calculated  values  for  the  highest  excited  state 
were  not  included  because  no  states  were  characterized  between  it  (at  3.96  eV) 
and  the  next  lowest  state  (at  1.74  eV).  As  mentioned  above,  the  highest  excited 
state  characterized  was  the  initial  converged  reference,  which  provided  an  initial 
wave  function  to  use  in  probing  for  lower-lying  states.  The  focus  of  further  char¬ 
acterization  was  to  survey  and  characterize  the  region  lying  between  0  and 
1.75  eV  of  the  ground  state  as  shown  in  Figure  2. 

As  can  be  seen  in  Figure  2,  the  calculated  values  of  the  electric  field  gradient 
at  the  central  atom  in  the  fragment  range  from  -0.0008  to  0.0153  a.u.  If  only  the 
ground  state  and  the  six  states  within  1.5  eV  of  the  ground  state  are  considered, 
the  range  is  considerably  smaller  (0.0091  to  0.0121  a.u.).  These  state-to-state  vari¬ 
ations  in  Bens  can  be  compared  to  state-to-state  variations  of  0.0159  to 
0.0211  a.u.,  -0.0165  to  0.0163  a.u.,  and  -0.0329  to  0.0526  a.u.  for  fragments  of 
69  [4],  81  [5],  and  87  [5]  atoms.  The  experimentally  observed  range  of  values  is 
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(2:0.0044-0.0049  a.u.  [22]).  Since  a  significant  trend  toward  the  experimentally 
observed  values  is  not  found  as  the  size  of  the  fragment  increases  from  69  to 
135  atoms  and  since  the  value  of  q  is  so  small  for  beryllium,  refinements  to  the 
present  treatment  may  be  required  for  highly  accurate  calculation  of  this  quan¬ 
tity.  Such  refinements  may  include  increasing  the  number  of  atoms,  including 
electron  correlation  effects,  extending  the  basis  set  on  the  central  atom,  employ¬ 
ing  an  all-electron  treatment  on  the  central  atom,  and  investigating  the  effects  of 
geometry  optimization  and  surface  states. 

An  additional  consideration  for  accurate  calculations  of  the  electric  field  gra¬ 
dient  as  discussed  in  an  earlier  study  of  a  cluster  of  63  beryllium  atoms  [3]  is  a 
correction  for  the  Sternheimer  effect  [1].  The  Sternhcimer  effect  takes  into  ac¬ 
count  polarization  of  the  core  (Is)  electrons  of  the  Be  atom.  This  effect  results  in 
a  16-18%  increase  in  magnitude  in  the  electric  field  gradient  which,  although  in 
the  wrong  direction  in  the  present  case,  does  not  qualitatively  alter  the  agree¬ 
ment  with  experiment. 

As  can  also  be  seen  in  Figure  2,  the  quadruple  moment  [23]  is  clearly  more 
sensitive  to  electron  configuration  since  the  principal  contributions  come  from 
the  outer  orbitals.  This  is  due  in  part  to  the  definition  of  the  quadrupole  moment 
as  a  difference  of  second  moments.  This  effect  can  be  seen  by  comparing  the 
state-to-state  variations  of  the  quadrupole  moment  in  the  figure  to  those  ob¬ 
served  for  the  nuclear-electrostatic  potential,  the  diamagnetic  shielding  constant, 
and  the  square  roots  of  second  moments  of  charge  discussed  above.  This  effect  is 
also  consistent  with  results  from  studies  on  smaller  fragments  [2-5]. 

Mulliken  population  analyses  [14]  have  been  carried  out  on  the  symmetry- 
distinct  groups  of  atoms  for  the  low-lying  states  reported  in  Table  I.  The  maxi¬ 
mum  state-to-state  variations  for  the  net  charges  were  found  to  be  less  than 
0.05  electrons  for  all  but  4  of  the  20  symmetry-distinct  groups  of  atoms  [4],  In  the 
remaining  4  groups  of  atoms,  the  maximum  state-to-state  variation  was 
0.1  electrons.  The  small  state-to-state  variations  are  consistent  with  the  earlier 
studies  on  bulk  fragments  ranging  from  51  to  87  atoms  [2, 4, 5]  and  are  indicative 
of  bulk  behavior  in  that  the  relative  populations  of  the  symmetry-distinct  sets  of 
atoms  are  essentially  independent  of  electronic  state. 

The  averages  of  net  charges  from  the  Mulliken  population  analyses  calculated 
for  the  low-lying  converged  states  are  shown  in  Figure  3  for  the  20  symmetry- 
distinct  groups  of  atoms  as  function  of  the  radii  of  the  coordination  shells  con¬ 
taining  the  respective  groups  [4],  As  can  be  seen  in  the  figure,  the  net  charge 
approaches  zero  at  the  center  of  the  cluster,  which  is  as  expected  in  the  neutral 
bulk  environment.  The  average  net  charges  on  the  cer  ral  Be  atom  for  bulk  frag¬ 
ments  of  51,  57,  69,  81,  87,  and  135  atoms  are  0.32,  0.40,  0.28,  0.15,  0.16,  and 
0.0  electrons,  respectively,  indicating  a  convergence  of  this  quantity  to  the  bulk 
value  in  a  fragment  of  135  atoms. 

Figure  3  also  shows  a  buildup  of  positive  charge  in  the  outer  region  of  the 
fragment  complemented  by  a  depletion  in  the  region  between  the  outer  and  cen¬ 
tral  regions.  This  may  be  interpreted  as  a  surface  effect  in  that  electrons  are  be¬ 
ing  drawn  from  the  surface  to  the  middle  region  of  the  fragment  so  as  to  be  in  a 
region  with  more  positive  charge.  It  may  also  be  an  indication  that  given  geomet¬ 
ric  flexibility,  the  surface  might  reconstruct  to  minimize  charge  separation. 
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In  addition  to  net  charges  on  the  symmetry-distinct  sets  of  atoms,  gross  atomic 
orbital  populations  were  also  computed.  On  the  central  Be  atom,  the  average 
values  of  the  s  and  p  gross  atomic  orbital  populations  were  found  to  be  0.56  and 
1.45  electrons,  respectively.  The  ratio  of  p:s  gross  atomic  orbital  populations  is 
consistent  with  the  valence  electron  density  analysis  obtained  from  X-ray  dif¬ 
fraction  studies  [24],  which  have  been  interpreted  to  indicate  s p*  hybrid  bonding 
in  beryllium  metal. 

The  molecular  valence  charge  densities  in  the  XY  (perpendicular  to  the  prin¬ 
cipal  Cj  axis  of  symmetry),  XZ,  and  YZ  planes  for  the  ground  state  of  the  frag¬ 
ment  are  shown  in  Figure  4.  The  plots  display  the  homogeneity  of  the  charge 
distribution  around  the  central  Be  atom  when  compared  to  charge  densities  cal¬ 
culated  for  fragments  of  81  and  87  atoms  [5],  thus  indicating  the  onset  of  a  con¬ 
stant  charge  distribution  around  the  central  Be  in  clusters  as  small  as  80  atoms. 

Evidence  of  p  participation  in  the  bonding  as  indicated  by  the  population 
analyses  and  X-ray  diffraction  studies  [24]  can  also  be  seen  in  Figure  4  by  exam¬ 
ining  the  anisotropy  of  the  bonding  in  the  XZ  and  YZ  planes  [Figs.  4(b)  and  4(c), 
respectively].  The  noncircular  shaped  contours  are  indicative  of  participation  of 
p  orbitals  in  the  interatomic  interactions. 

Evidence  of  electronic  charge  density  maxima  in  regions  between  nuclei  can 
be  seen  in  Figure  4(a)  in  the  XY  plane.  These  maxima,  however,  are  not  in  the 
proper  location  in  the  fragment  to  correspond  to  the  electron  density  maxima 
predicted  to  be  located  near  tetrahedral  hole  sites  in  the  hexagonal  closed-pack 
lattice  by  electron  density  deformation  and  valence  charge  density  plots  from  an 
X-ray  diffraction  study  [24]  and  a  Hartree-Fock  band  theory  study  [7].  The  cur¬ 
rent  studies  are  consistent  with  local  density  functional  band  theory  (ldfb)  stud¬ 
ies  [9,10],  which  also  do  not  predict  density  maxima  at  the  tetrahedral  hole  sites. 

The  density  of  states  (dos)  for  the  ground  state  of  the  fragment  can  be  seen  in 
Figure  5.  We  find  that  the  Dos  for  the  canonical  Hartree-Fock  orbitals  is  zero  at 
an  energy  immediately  following  the  Fermi  level.  This  is  consistent  with  an  ear¬ 
lier  study  by  Monkhorst  showing  that  the  dos  at  the  Fermi  level  is  zero  in  canoni¬ 
cal  crystal  Hartree-Fock  studies  on  partially  filled  band  systems  [25].  It  is  also 
consistent  with  a  hfdos  study  by  Cox  and  Bauschlicher  on  Ah»  in  which  a  gap  in 
the  dos  appeared  near  the  Fermi  level  [26].  In  the  present  study  the  dos  using 
canonical  virtual  orbitals  is  zero  just  above  the  Fermi  level  since  the  Fermi  level 
is  taken  as  the  energy  of  the  highest  occupied  state.  In  contrast,  the  dos  employ¬ 
ing  the  tvo  orbitals  is  nonzero  in  the  Fermi  level  region,  consistent  with  the  fact 
that  beryllium  is  a  conductor. 

The  general  features  observed  in  the  ivo  Dos  are  a  dip  near  the  Fermi  level,  a 
maximum  in  the  occupied  region  between  -0.25  and  -0.30  a.u.,  and  a  rapid  rise 
in  the  unoccupied  region  leading  to  a  broad  maximum  starting  at  0.35  a.u.  These 
general  features  indicate  a  rapid  decrease  of  states  in  the  energy  region  before 
the  Fermi  level  followed  by  a  rapid  increase  in  unoccupied  states.  These  features 
were  found  consistently  in  dos  generated  using  a  range  of  energy  increments 
from  0.035  to  0.080  a.u. 

The  density  of  states  is  not  a  direct  experimentally  observable  quantity,  but  in¬ 
terpretation  of  photo-electron  and  X-ray  emission  and  absorption  spectra  can 
provide  insight  into  the  dos  [27],  Detailed  comparisons  between  the  ivodos  and 


experiment  must  take  into  account  other  factors  such  as  intensity  matrix  ele¬ 
ments,  final-state  relaxation  effects,  and  electron  correlation.  Notwithstanding 
those  considerations,  X-ray  emission  and  absorption  studies  in  the  present  case 
show  evidence  of  a  dip  near  the  Fermi  level  [28],  a  rapid  rise  in  unoccupied 
states  [28],  and  fine  structure  in  the  occupied  dos  toward  the  low  end  of  the  en¬ 
ergy  range  [29].  These  features  are  consistent  with  the  current  tvo  dos  and  previ¬ 
ous  dos  analyses  calculated  for  beryllium  by  Hartree-Fock  band  (hfb)  [7],  local 
density  functional  band  (ldfb)  [8,9],  and  orthogonalized  plane  wave  (opw)  [10] 
methods.  The  dip  near  the  Fermi  level  is  also  consistent  with  an  additional  ear¬ 
lier  opw  study  by  Herring  and  Hill  [11]. 

As  ca  be  seen  in  Figure  5,  the  bandwidth  of  the  occupied  region  of  the  dos  is 
0.61  ax.  This  can  be  compared  to  a  bandwidth  of  0.7  a.u.  estimated  from  the 
dos  in  the  hfb  study  [7].  The  experimental  bandwidth  as  estimated  from  experi¬ 
mental  X-ray  emission  spectra  is  0.46  a.u.  [30].  The  larger  theoretical  bandwidth 
is  consistent  with  earlier  Hartree-Fock  studies  and  has  been  shown  to  decrease 
when  electron  correlation  effects  are  included  [25, 26].  Earlier  estimates  of  band- 
widths  range  from  0.39  to  0.45  a.u.  depending  on  whether  the  ldfb  [8,9]  or  opw 
[10,11]  model  is  used. 


Conclusions 

An  ab  initio  study  on  a  bulk  fragment  of  135  atoms  has  been  carried  out  to 
investigate  the  electronic  structure  and  density  of  states  in  beryllium  metal.  Con¬ 
vergence  to  the  bulk  limit  occurs  for  the  net  population  on  the  central  Be  atom. 
The  onset  of  convergence  to  the  bulk  limit  is  observed  for  the  binding  energy  per 
atom.  In  addition,  ionization  potentials  for  the  lowest  lying  states  arc  within  10% 
of  the  bulk  experimental  work  function.  The  density  of  states  employing  im¬ 
proved  virtual  orbitals  compares  well  with  experimental  evidence.  First-principle 
studies  on  fragments  of  this  size  are  shown  to  give  valuable  insight  into  the  de¬ 
scription  of  the  electronic  structure  of  bulk  systems  and  are  in  agreement  with 
‘■'•errintions  obtained  both  from  experiment  and  from  band  theory  models. 
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Abstract 

Ab  initio  spin-orbit  full  configuration  interaction  calculations  in  the  context  of 
relativistic  effective  core  potentials  are  reported  for  the  weakly  bound  metal  dimer 
LiBe,  a  3-valcncc-elcctron  system.  The  effects  of  basis  set  on  the  energies  of  valence 
and  Rydberg  states  of  the  cluster  arc  discussed,  as  arc  the  effects  of  configuration 
space  selection  on  the  energy  of  the  latter  states.  Results  at  the  dissociative  limit  arc 
compared  to  the  experimental  atomic  spectra.  Potential  energy  curves  and 
spectroscopic  constants  are  presented  for  the  ground  state  and  fourteen  excited 
states,  which  includes  the  Li  and  Be  2p  valence  states,  the  Li  3s,  3p,  3d  and  4s  Rydberg 
states,  as  well  as  three  low-lying  states  of  the  molecular  cation. 


37 


I.  Introduction 

There  is  great  current  interest  in  the  structure  and  chemical  properties  of 

small  metal  clusters.  The  reasons  for  the  increased  research  activity  in  this  area  arc 

several-fold.  First,  recent  advances  in  the  generation  of  these  species,  particularly  by 
laser  vaporization  techniques,  have  made  it  feasible  to  study  them  experimentally. 3  At 
the  same  time,  the  development  of  new  theoretical  methods,  combined  with  the  latest 
advances  in  computing  technology,  have  made  possible  the  treatment  of  all  types  of 
metal  clusters,  including  those  containing  heavy  elements,  at  a  post-Hartree-Fock  level 

of  theory. 4-6  Such  parallel  advances  in  experiment  and  theory  provide  an  ideal  setting 
for  characterizing  these  systems. 

An  accurate  theoretical  study,  however,  requires  the  availability  of 

computational  methods  which  are  sophisticated  enough  to  generate  reliable  electronic 
spectra,  potential  energy  curves  and  wave  functions  while  still  remaining  tractable.  One 
of  the  most  powerful  general  methods  to  date  is  full  configuration  interaction  (full  Cl).  It 
is,  however,  computationally  demanding  and,  depending  on  the  number  of  electrons, 
impossible  to  perform  on  heavy  elemental  clusters  because  of  the  large  number  of  terms 
which  must  be  stored  and  manipulated.  The  method  employed  in  this  study  replaces  the 

core  electrons  of  each  atom  with  averaged  relativistic  effective  potentials  (AREP's), ^  thus 
reducing  the  electrons  in  the  valence  region  to  a  number  small  enough  such  that  a  full 
Cl  treatment  is  feasible.  In  this  case,  the  only  errors  introduced  arise  from  the  AREP's 
and  the  finite  basis  set.  If  the  AREP's  used  arc  determined  to  be  accurate,  then  the  error 
is  attributable  almost  solely  to  inadequacies  in  the  basis  set.  It  is  then  possible  to 
determine  the  most  appropratc  basis  set  for  the  accurate  calculation  of  the  electronic 
spectra  of  the  system. 

This  paper  reports  the  effects  of  basis  set  changes  on  the  electronic  energies 
of  the  valence  and  Rydberg  states  of  LiBc  using  a  full  Cl  treatment  that  incorporates  the 
spin-orbit  operator  explicitly  in  the  Hamiltonian  (SOCI).  The  incorporation  of  a  spin- 
orbit  operator  increases  the  accuracy  of  the  calculated  results.7  The  core  electrons  of 
both  Li  and  Be  are  represented  by  AREP's,  thereby  reducing  the  number  of  electrons  in 
the  valence  region  to  exactly  three.  Potential  energy  curves  and  spectroscopic  constants 
for  the  ground  and  several  excited  stales  of  LiBc  and  LiBc+  arc  presented,  and,  where 
possible,  the  results  arc  compared  with  experiment. 
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II.  Calculations 

The  SOCI  procedure  involves  a  Hamiltonian  matrix  that  is  constructed  in  terms 
of  elements  involving,  in  addition  to  the  terms  in  the  nonrelativistic  Hamiltonian 
operator  for  the  valence  electrons,  averaged  relativistic  effective  potentials  (AREP),  a 
spin-orbit  hamiltonian  Hso  and  determinanta!  wavefunctions  for  valence  electrons 
defined  through  standard  L-S  coupling  SCF  procedures.^  The  REP's  are  defined  using 
shape-consistent  nodelcss  pseudo-orbitals*®  extracted  from  numerical  two-component 
Dirac-Fock  atomic  spinors*  *  using  the  procedure  of  Lee  et  al.*2  The  pseudo-orbitals  arc 
derived  such  that  nodes  in  the  respective  atomic  spinors  are  removed  by  requiring  that 
up  to  five  derivatives  be  matched  at  a  distance  from  the  atomic  nucleus  beyond  which  the 
nodeless  orbital  is  exactly  the  same  as  the  true  Dirac-Fock  solution.*®  Such  pseudo¬ 
orbitals  are  then  used  to  define  individual  REPs  in  j-j  coupling.  The  j  ±  ~  components  arc 

then  used  to  form  weighted  average  REP  (AREP)  operators  and  spin-orbit  (SO) 
operators. Because  the  SOCI  calculations  result  in  wavefunctions  that  transform 
according  to  one  of  the  irreducible  representations  of  the  C2v'  double  group,  and  since 
LiBc  possesses  an  odd  number  of  electrons,  all  possible  states  are  extracted  from  a  single 
Hamiltonian  matrix  of  E  symmetry.^  This  is  a  clear  advantage  in  the  full  Cl  procedure  in 
that  energies  and  wavefunctions  of  all  electronic  states  for  a  given  intemuclcar 
separation  can  be  extracted  as  successive  eigensolutions  of  one  matrix  provided  the 
atomic  basis  sets  are  adequate. 

The  AREP's  used  in  this  work  describe  the  Is  electrons  of  both  Li  and  Be.*  3 
Basis  sets  of  contracted  Gaussian  functions  were  optimized  at  the  atomic  level  and  used  to 
describe  the  valence  states  of  both  atoms.  These  sets  are  given  in  Table  I.  In  the  case  of 
Li,  one  diffuse  Gaussian  function  was  included  initially  in  the  basis  to  describe  the 
Rydberg  3s  and  4s  states.  The  basis  functions  for  the  Li  3p  and  3d  Rydberg  states,  as  well 
as  improved  basis  functions  for  the  Li  3s  and  4s  Rydberg  states,  were  composed  of 
Gaussian  expansions  of  Slater-type  orbitals  (STO-6G).*^  The  exponents  of  the  latter 
functions  were  scaled  to  describe  the  corresponding  orbitals  of  the  Li  atom.  The  scaling 
parameter  was  determined1  ^  from  the  quantum  defect  using  experimental  ionization 
potentials  and  the  corresponding  atomic  energies  of  the  Rydberg  states  of  the  Li  atom. 

For  each  Rydberg  orbital,  £  =  1/n*  ,  where  £  is  the  exponent  for  the  Slater  orbital 
representing  a  given  Rydberg  state,  n*  =  n  -  5  and  S  is  the  quantum  defect  of  the  orbital 
of  principal  quantum  number  n.  Here,  n*  =  [  R/(AEj0n  -AEr)]*/^  and  AEj0n  and  AEr 
arc,  respectively,  the  experimental  atomic  ionization  and  excitation  energies 
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corresponding  to  a  particular  Rydberg  state.16  R  is  the  Rydberg  constant.  (R  =  109734 
cm"').  Given  the  experimental  energies,  it  is  possible  to  calculate  £  for  each  Rydberg 

orbital.  Exponents  a  for  such  values  of  £  are  obtained  from  those  for  £=1  by  multiplying 

each  a  in  the  Gaussian  expansion  of  Slater-type  orbitals  by  £2.14  The  exponents  and 
coefficients  of  the  resulting  Rydberg  STO-6G  functions  are  given  in  Table  II. 

Energies  corresponding  to  intemuclear  separations  from  3.5  to  10  bohr  were 
calculated  for  the  ground,  excited  valence  and  ionic  states.  The  calculated  energies  of  the 
Rydberg  states  correspond  to  intemuclear  distances  from  3.5  to  7  bohr;  specifically,  3.5, 
3.95,  4.5,  4.896,  5.5,  5.85,  6.5,  and  7.0.  The  energies  corresponding  to  distances  greater 
than  7  and  up  to  10  bohr  for  these  states  were  derived  from  a  fit  to  a  Hulburt-Hirshfelder 
(H-H)  analytical  form.17  The  energies  of  all  the  states  of  the  neutral  and  cationic  systems 
were  also  calculated  for  an  intemuclear  distance  of  40  bohr.  In  addition,  calcuations 
were  completed  for  R=106  bohr  for  the  cation.  Spectroscopic  constants  for  the  ground, 
excited  valence,  Rydberg  and  ionic  states  were  calculated  to  second  order  of  perturbation 
theory  by  fitting  the  energies  of  each  to  a  high-degree  (6th  or  higher)  Taylor  series. 

In  certain  cases  a  selection  scheme  was  used  to  obtain  the  spatial 
configurations  comprised  of  double-group  symmetry-adapted  MO's  that  were  retained  in 
the  SOCI  wavefunction.  Two  cases  were  treated  to  determine  the  effects  of  such 
configuration  selection  on  the  energies  of  the  Rydberg  and  valence  states  of  the  system. 
The  first  involved  inclusion  of  the  complete  set  of  valence  configurations  plus  all 
Rydberg  configurations.  The  second  case  included  all  valence  and  Rydberg 
configurations  as  well  as  all  configurations  corresponding  to  single  excitations  relative 
to  the  Rydberg  configurations.  In  each  instance,  the  Cl  wavefunction  was  calculated 
using  configurations  derived  from  the  la22a  SCF  MO's  of  the  LiBe  ground  state  plus 

orthonormal  virtual  (unoccupied)  MO's,  and  the  results  were  compared  to  the 
corresponding  full  SOCI  values. 

The  dissociative  limit  was  assumed  to  be  40  bohr,  at  which  point  the  resulting 
theoretical  values,  were  compared  to  atomic  spectra.16.  All  calculations  were  performed 
on  a  Cray  Y-MP  supercomputer  using  codes  that  are  part  of  the  COLUMBUS  suite  of 
programs. 9’  18  The  calculations  involving  the  largest  basis  set  (56  basis  functions) 
required  approximately  25  hours  of  epu  time  to  construct  and  extract  as  many  as  50  roots 
of  the  Hamiltonian  matrix  comprised  of  30,800  double  group  symmetry  adapted 
configurations  and  113,960  determinantal  wave  functions.  The  Hamiltonian  matrix 
required  approximately  0.75  gigabytes  of  disk  storage. 
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III.  Results  and  Discussion 
A.  Electronic  Spectrum 

The  theoretical  spectrum  of  LiBe  at  dissociation,  including  the  3s,  3p,  3d  and  4s 
Rydberg  states  of  Li,  each  represented  by  an  STO-6G  function  (basis  set  6  in  Table  I),  is 
given  in  Table  III.  State  1  of  Table  III  is  the  ground  state,  states  2  through  4  correspond  to 
the  Li  2p  (2P)  states,  while  states  5  through  13  correspond  to  the  Be  2p  (3P)  states.  The  Li 
Rydberg  3s,  3p,  3d  and  4s  states  arc  14,  15  through  17,  18  through  22  and  state  23, 
respectively.  Although  not  shown  in  Table  III  there  are  a  total  of  27  additional  states  that 
correspond  to  double  excitation  states  of  the  LiBe  system.  Simultaneous  single  excitations 
of  both  Li  and  Be  electrons  give  rise  to  these  states,  which  lie  at  approximately  4.61  eV 
and  have  energies  that  are  sums  of  the  single  excitation  energies  of  Li  and  Be.  Also 
shown  are  low-lying  states  of  LiBe+.  The  1st  and  ll1*1  of  these  are  the  ionization 
potentials  of  Li  and  Be,  respectively,  while  the  2nd  through  10th  correspond  to  the 
ionization  of  Li  plus  the  Be  2s2p  (3P)  excitation.  All  of  the  calculated  values  are  in 
excellent  agreement  with  experiment,  the  largest  error  being  0.05  eV  for  both  the  3d 
states  and  ionization  potential  of  Li. 

The  potential  energy  curves  for  LiBe  were  calculated  using  basis  set  6.  Figure 
1  shows  the  curves  calculated  for  the  ground,  valence  and  cation  states  and  their 
corresponding  atomic  states  at  dissociation.  The  three  highest-lying  curves  shown 
correspond  to  the  three  lowest  states  of  LiBe+-  The  dissociative  limit  of  the  first  of  these  is 
Li+  +  Be.  The  second  and  third  correspond  to  dissociation  to  Li+  +  Be(3P).  The  states  lying 
between  -1.13  and  -1.06  hartrees  are  shown  on  an  expanded  scale  in  Figure  2. 

Calculated  atomic  spin-orbit  splittings  arc  in  excellent  agreement  with 
experimental  values16  for  both  the  neutral  and  ionic  systems.  The  greatest  discrepancy 
is  0.1  meV  (Table  III).  The  largest  spin-orbit  splitting  calculated  is  0.2  meV  between  both 
the  neutral  2S  +  3Pi  and  2S  +  3P2  and  the  ionic  *S  +  3Pj  and  JS  +  3P2  states.  The 
experimental  splittings  are  0.3  meV  for  each.  The  rest  of  the  states  have  splittings  in  the 
range  0.0-0. 1  meV  as  both  calculated  and  observed. 

The  spin-orbit  splittings  corresponding  to  the  equilibrium  internuclear 
separations  are  also  small.  The  A2!!,  C2n,  b4Z',  2n(3p),  2A(3d)  and  2n(3d)  states  of  neutral 
LiBe  are  each  split  into  two  spin-orbit  components.  ‘  The  splittings  are  0.08,  0.0,  0.0,  0.008, 
0.0  and  0.03  meV,  respectively.  The  a4n  state  contains  four  spin-orbit  components.  Each 
is  separated  by  0.04  meV  from  the  state  or  states  lying  closest  to  ii  in  energy,  making  the 
total  spin-orbit  energy  range  0.12  meV.  The  B3fl  state  of  LiBe+  is  split  into  three  spin- 
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orbit  states.  There  is  a  splitting  of  0.22  meV  between  the  lowest  component  and  the 
higher  two,  which  are  degenerate. 

The  equilibrium  intemuclear  separation  (Re),  harmonic  frequency  (oje). 
anharmonicity  constant  (o)exe),  rotational  constant  (Be),  vibration-rotation  coupling 
constant  (ae),  zero-point  vibrational  energy  (Go),  and  minimum-to-minimum  electronic 
excitation  energy  (Te)  for  the  ground,  excited  valence,  Rydberg  and  ionic  states  appear  in 
Table  IV.  Spectroscopic  constants  for  the  ground  state  are  in  excellent  agreement  with 
those  calculated  previously  using  a  numerical  solution  to  the  nuclear  motion  Schroedinger 
equation. Most  of  the  states  have  similar  wc  values  with  the  exceptions  of  the  C2n,  D2L+, 
A2n  and  b4X*  states  The  first  two  have  a  low  harmonic  frequency  due  to  their  broad, 
shallow  wells.  The  latter  two  have  wells  that  are  relatively  deep  -  which  accounts  for  their 
higher  frequencies  -  and  similar  in  shape  and  thus  have  very  close  values  of  ue.  Most  of 
the  rotational  constants  have  comparable  values  (between  0.5  and  0.7  cm*1),  as  expected, 
since  the  equilibrium  intemuclear  separations,  Re  values,  lie  between  5  and  6  bohr.  One 
exception  is  the  D2I+  state,  whose  Be  value  is  low  (0.33  cm'1)  compared  to  the  others  due  to 
its  relatively  long  equilibrium  intemuclear  separation  (6.73  bohr).  The  other  two 
exceptions  are  the  A2n  and  b4X~  states  whose  Revalues  are  closer  to  4  bohr,  resulting  in 
larger  Be  values.  The  equilibrium  rotation-vibration  coupling  constants  are  also  similar 
for  all  states  except  the  D2X+  and  2n(3d),  which  have  negative  values  of  ae.  The  latter 
curve  actually  starts  out  at  small  R  as  a  repulsive  state.  However,  the  avoided  crossing 
between  it  and  the  2n(3p)  state  causes  the  curve  to  swing  up  and  thus  become  more 
harmonic.  The  low  anharmonicity  of  the  state  results  in  a  negative  ae.  The  former  state 
also  has  a  negative  value  of  ae  In  this  case,  the  D2Z+  state  would  have  had  a  much  deeper 
well.  However,  the  avoided  crossing  between  it  and  the  2X+(3s)  Rydberg  state  causes  the 
former  curve  to  swing  down,  resulting  in  a  long,  shallow  minimum  with  a  low 
anharmonicity.  Again,  all  of  the  states  have  comparable  Go  values  with  the  exception  of 
the  C2n  and  D2X+  states,  which  have  low  Go  (96  and  85  cm*1,  respectively).  This  is  due  to 
the  fact  that  these  two  states  are  very  broad  and  shallow.  Te  values  show  that  there  are  two 
low-lying  excited  valence  states  (A2n  and  B2X+).  The  rest  of  the  states  are  more  than  10.000 
cm*1  higher  in  energy  relative  to  the  ground  state.  It  is  noted  that  the  b4X*  valence  state 
is  close  in  energy  to  the  Rydberg  2L+(3s)  state,  the  energy  difference  being  approximately 
1000  cm*1.  It  is  therefore  possible  that  the  Rydberg  and  higher  valence  states  interact 
This  is,  in  fact,  the  case  between  theD2Z+and  2Z+(3s)  states,  which  undergo  an  avoided 
crossing.  Therefore,  an  accurate  spectral  study  requires  the  inclusion  of  Rydberg 
functions  in  the  basis  set. 

Table  IV  also  lists  the  dissociation  energies  (De)  for  all  of  the  states  calculated. 


42 


The  ground  state  of  LiBe  is  bound  by  0.25  eV.  This  is  in  agreement  with  previous 
calculations.  18-21  The  latter  reference  further  points  out  that  the  bonding  is  largely  due 
to  van  der  Waals  interactions  and  substantiates  this  by  showing  that  the  dissociation 
energy  of  LiBe  is  proportional  to  l/rc^.  The  remaining  states  of  both  the  neutral  and 
ionic  systems  are  weakly  bound,  with  the  exceptions  of  the  A2n  and  b4Z*  valence  states, 
which  are  bound  by  1.50  eV  and  2.03  eV,  respectively,  and  the  valence  4£+  (I)  state.  The 
latter  involves  a  repulsive  interaction  (Figure  2).  Despite  weak  bonding  in  most  cases,  all 
of  the  excited  states  for  the  neutral  system,  as  well  as  the  two  lower-lying  states  of  the 
ion,  exhibit  more  bonding  than  the  ground  state.  The  equilibrium  intemuclear  distances 
(Re's)  of  all  the  states  spread  over  a  fairly  large  range,  approximately  4.0  -  7.0  bohr. 

As  shown  in  Figure  1,  the  first  three  states  calculated  for  the  neutral  system 
and  all  three  states  of  the  ion  are  simple,  i.e.,  there  are  relatively  few  interactions  among 
these  due  to  the  large  energy  differences  involved.  Where  the  curves  approach  each 
other,  only  two  crossings  -  between  the  X  and  A  states  at  3.6  bohr  and  between  the  A  and  B 
states  of  LiBe  at  6.8  bohr  -  are  calculated.  The  major  area  of  interaction  occurs  among  the 
upper  valence  and  Rydberg  states.  This  area  is  shown  on  a  larger  scale  in  Figure  2. 
Previous  calculations18  agree  very  well  with  those  reported  here  for  the  valence  states. 
However,  that  work  focussed  on  the  valence  states  of  LiBe  and  therefore  the  basis  set  did 
not  include  Rydberg  functions.  Because  of  the  interactions  among  the  upper  valence 

and  the  Rydberg  states,  an  accurate  representation  of  the  valence  states  lying  higher  in 
energy  than  the  a4n  state  could  not  be  made  without  inclusion  of  the  Rydberg  functions 
to  the  basis  set.  Discrepancies  between  the  upper  valence  curves  presented  here  and 

those  of  the  previous  work  are  a  consequence  of  the  differences  in  the  basis  sets  used. 

As  stated  previously,  the  most  interesting  region  in  Figure  1  appears  on  an 

expanded  scale  in  Figure  2.  There  are  three  avoided  crossings  shown  here.  The  first 

occurs  at  approximately  4.8  bohr  and  is  between  the  D2L+  and  the  Rydberg  2I+  (3s)  states. 
This  is  the  reason  for  the  large  anharmonicity  seen  in  the  latter  (coexe  is  9.2  cm*1)  The 
second  avoided  crossing  occurs  at  7.3  bohr  and  is  between  the  2I+  (3d)  and  2I+  (4s)  states, 
leading  to  a  higher  anharmonicity  for  these,  also.  However,  since  the  avoided 
interaction  is  farther  from  the  equilibrium  intemuclear  separation  for  both,  their  cocxe 
values  (5.0  and  6.0  cm*1,  respectively)  are  smaller  relative  to  those  of  the  state  involved 
in  the  previous  avoided  crossing.  The  third  avoided  interaction  occurs  at  approximately 
5.8  bohr  between  states  2n  (3p)  and  2FI  (3d).  Again,  toexe  is  large  (6.7 

cm*1)  for  the  former  state  due  to  this  interaction.  If  the  crossing  were  allowed,  the  latter 

state  would  be  repulsive.  However,  as  a  result  of  the  avoided  interaction  with  the  2n  (3p) 

state,  it  curves  upward,  forming  a  shallow  well  with  a  relatively  small  coexc  (3.1  cm*1). 
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It  was  mentioned  above  that  the  region  greater  than  7  bohr  in  Figure  1  and 
corresponding  to  the  Rydberg  states  contains  energies  which  were  approximated  using  a 
H-H  analytical  form.  This  function  is  excellent  for  a  general  characterization  of  each 
curve,  especially  in  the  region  between  the  equilibrium  intemuclear  distance  and  the 
dissociative  limit.  However,  the  description  is  expected  to  be  less  accurate  in  cases  where 
the  interactions  among  curves  results  in  a  deviation  from  the  usual  curve  behavior,  as  in 
the  case  of  avoided  crossings  For  most  of  the  curves  in  Figure  1  this  approximation  holds 
very  well.  For  example,  the  difference  between  the  H-H  energy  and  the  calculated 
energy  of  the  ground  state  at  8  bohr  is  0.2  mhartrees  (-1.2057  vs  -1.2059,  respectively). 

At  10  bohr  it  is  0.3  mhartrees  (-1.2046  vs  -1.2049).  In  the  region  corresponding  to 
intemuclear  separations  between  7  and  10  bohr  there  are  five  cases  where  the  curves 
either  cross  or  avoid,  and  thus  five  possible  areas  where  the  H-H  approximation  may 
break  down.  The  lowest-energy  of  these  is  the  crossing  between  the  a4n  and  the 
repulsive  4£+(I)  states  at  approximately  9  bohr.  This  crossing  is  predicted  by  the  H-H 
approximation  and  verified  by  calculation.  The  agreement  between  the  two  methods 
results  from  the  fact  that  even  though  the  curves  cross,  the  crossing  does  not  alter  their 
shape.  The  second  case  involves  the  2Z+(3p)  and  2n(3p)  states.  These  curves  are 
calculated  to  cross  at  about  7  bohr.  After  the  crossing,  the  latter  swings  down  but  must 
curve  up  again  since- it  cannot  dissociate  to  the  3s  limit.  It  is  then  possible  that  it  may 
cross  the  2I+  state  again  or  parallel  it  toward  dissociation.  The  second  choice  is  more 
probable.  This  curve  cannot  dissociate  to  a  3d  asymptote,  and  the  likelihood  that  it  will 
cross  the  2I+  state  again  decreases  as  the  dissociative  limit  is  approached.  The  third  case 
is  the  crossing  between  the  2I7(3d)  and  2A  (3d>  states  at  about  7.5  bohr.  The  avoided 
crossing  between  the  2n(3d)  and  the  2n(3p)  states  at  approximately  6  bohr  discussed 
earlier  causes  the  former  curve  to  swing  upward,  crossing  the  2  A  (3d)  state  again.  The 
curve  continues  upward  to  cross  the  2I+(3d)  state  (the  fourth  case)  before  dissociating. 
The  latter  crossing  occurs  because  the  calculated  avoided  interaction  between  the  2I+(3d) 
and  2£+(4s)  states  discussed  previously  causes  the  former  state  to  dip  down  before 
dissociating.  The  last  area  of  interest  is  the  avoided  crossing  between  the  2£+(3d)  and 
2I+(4s)  states.  Although  the  curves  actually  avoid  at  about  7.5  bohr,  this  behavior  is 
easily  predicted  at  approximately  6,5-7  bohr  from  the  slope  of  the  curves  and  the  fact  that 
a  crossing  between  them  is  forbidden  due  to  symmetry.  The  behavior  of  the  higher- 
lying  curves  is  expected  to  be  approximated  very  accurately  by  the  H-H  function  because 
the  energy  differences  among  the  curves  is  large  enough  to  ensure  no  interactions 
between  7  bohr  and  dissociation. 

This  study  involved  the  calculation  of  25  states  of  the  LiBe  and  LiBc+  systems. 
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18  of  which  appear  in  Figure  1.  The  rest  are  spin-orbit  components  of  the  those  pictured 
and  are  not  resolved  on  the  energy  scale  shown.  To  date  no  studies  of  this  magnitude 
have  been  published  for  LiBe.  However,  an  experimental  and  a  theoretical  spectrum  has 
been  published  for  LiMg.22  Since  Mg  is  a  relatively  light  clement  belonging  to  group 
Ila,  this  system  is  similar  to  LiBe  and  therefore  the  spectrum  of  one  may  be  compared  to 
that  of  the  other.  The  first  five  states  calculated  for  LiMg  -  X2Z+,  A2n,  B2X+,  C2n  and  D2I+ 
-  are  very  similar  in  shape  and  relative  positions  to  the  corresponding  states  of  LiBe.  Bc 
values  for  the  last  two  LiMg  states  listed  are  0.4737  and  0.3583  cm'1,  respectively, 
compared  to  the  corresponding  values  of  0.459  and  0.337  cm'1  for  LiBe.  As  in  LiBe,  the 
X2X+,  C2n  and  D2E+  states  are  more  shallow  than  the  A2n  and  B2E+.  However,  our 
calculations  indicate  that  the  D2I+state  avoids  the  higher-lying  2L+(3s)  Rydberg  state, 
causing  a  slight  shoulder  in  the  former  at  approximately  5  bohr.  Although  it  is  unclear 
whether  a  similar  interaction  occurs  in  LiMg  because  the  atomic  Mg  3P  and 
excitations  arc  lower  in  energy  than  the  corresponding  ones  in  Be,  it  is  possible  that 
interactions  do  occur  among  the  lower-lying  Rydberg  states  and  the  higher-lying 
valence  ones  and  that  these  interactions  would  alter  the  shape  of  the  reported  curves.2  2 


B.  Basis  Set  Effects 

The  effects  of  basis  set  on  the  Li  Rydberg  3s  and  4s  states  were  studied  using 
five  different  sets,  each  of  which  varied  in  both  the  number  and  type  of  valence  and 
Rydberg  functions  comprising  them.  The  energy  levels  computed  with  these  basis  sets 
appear  in  Table  V.  All  values  are  in  eV  and  represent  energy  differences  between  the 

ground  and  excited  states.  The  basis  sets  are  described  in  Table  I. 

In  all  cases,  state  1  is  the  ground  state  and  states  2-4  correspond  to  Li  2p  (2P) 

excitations.  The  three  nearly  degenerate  roots  observed  are  a  result  of  the  incorporation 
of  spin-orbit  effects  into  the  calculation.  States  5-13  correspond  to  Be  2p  (3P)  excitations. 
Again,  the  presence  of  the  spin-orbit  operator  causes  the  appearance  of  nine  nearly 
degenerate  roots.  State  14  of  basis  sets  2,  3,  4  and  5  corresponds  to  the  Li  3s  Rydberg  state, 
while  state  15  of  basis  set  3  and  state  23  of  basis  set  5  correspond  to  the  Li  4s  Rydberg 
state.  The  higher-lying  states,  all  occurring  at  aproximately  4.39  eV  for  the  first  four 
basis  sets  mentioned,  correspond  to  poor  representations  of  higher-lying  Rydberg  states 
of  LiBe.  The  energies  are  slightly  improved  (from  4.389  to  4.392  eV)  as  the  valence  basis 
set  on  Li  is  increased  from  double  to  triple  zeta  and  one  Li  s  Rydberg  function  is  added  to 
the  basis  set.  This  clearly  indicates  that  basis  sets  1  through  4  provide  an  inaccurate 
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description  of  the  LiBe  electronic  spectrum  at  dissociation.  This  inadequate  description 
of  these  states,  either  by  a  lack  of  Li  Rydberg  functions  (basis  set  1)  or  an  incomplete  set 
of  Rydberg  functions  (basis  sets  2  through  4,  which  are  missing  functions  that  describe 
the  3p  and  3d  Rydberg  states)  causes  their  energies  to  appear  much  higher  than  they 
actually  are.  Because  the  SOCI  calculation  is  variational,  the  use  of  a  small  basis  set 
results  in  higher-than-actual  energy  values.  As  the  basis  set  is  improved,  the  energy  of 
the  states  would  decrease,  thereby  approaching  their  proper  values.  The  appearance  of 
the  Li  (3s)  Rydberg  state  at  3.8779  eV  (basis  set  2)  does  not  affect  the  energies  of  the 
remaining  states.  This  also  applies  to  the  appearance  of  both  the  Li  2  s  (3s)  and  2S  (4s) 
Rydberg  states  at  3.4179  and  4.3169  cV,  respectively  (basis  set  3).  The  energies  of  stales  14 
through  16  in  column  5  (basis  set  4)  arc  not  affected  relative  to  the  energies  of  the  same 
states  in  column  4  (basis  set  3)  despite  energy  differences  of  0.05  and  0.03  eV  between  the 
Li  3s  and  4s  Rydberg  states,  respectively,  in  these  two  columns.  In  other  words,  basis  set  4 
leads  to  a  lowering  of  0.05  eV  in  the  energy  of  state  14  relative  to  the  energy  of  this  state 
in  column  4,  which  corresponds  to  basis  set  3.  However,  the  energies  of  states  16  through 
18  are  only  0.00005  cV  higher  relative  to  the  energies  of  these  same  states  in  column  3. 
The  lowering  of  the  Rydberg  3s  and  4s  states  of  Li  to  their  proper  values  does  not  shift  the 
energies  of  the  remaining  Rydberg  states  to  their  correct  values  of  approximately  3.8  eV. 
The  error  in  these  energies  therefore  predominantly  arises  from  the  absence  of  Rydberg 
3p  and  3d  functions  in  the  Li  basis  set.  These  functions  would  give  rise  to  the  Li  2p  (3p) 
and  (3d)  states,  which  occur  at  3.8  eV.^  This  is,  in  fact,  the  case  for  basis  set  5,  which 
is  the  most  accurate  of  the  basis  sets  appearing  in  Table  V.  This  basis  set  contains  one 
STO-6G  function  to  represent  each  3s,  3p,  3d  and  4s  Rydberg  state  of  Li  (see  Table  I).  The 
resulting  energies  are  in  very  close  agreement  to  those  corresponding  to  basis  set  6  in 
Table  III. 

Despite  the  fact  that  basis  set  1  is  small  (double  zeta  quality  for  the  description 
of  both  s  and  p  valence  states  of  Li  and  Be),  the  energies  of  the  Li  2p  and  Be  3p  valence 
states  are  very  close  to  the  experimental  atomic  values,  the  energy  differences  being 
only  0.01  and  0.05  eV,  respectively.  Expansion  of  the  Li  valence  basis  set  from  double  zcia 
s  (basis  sets  1  and  2)  to  triple  zeta  s  (basis  sets  3  and  4)  lowers  the  energy  difference 
between  the  experimental  and  theoretical  values  corresponding  to  the  Li  2p  valence 
states  by  0.003  eV.  As  expected,  the  energies  of  states  5  through  13  are  not  affected  by 
this  change  since  they  correspond  to  Be  valence  states.  The  energies  of  the  Li  Rydberg 
states,  on  the  ether  hand,  are  too  high,  a  result  of  absent  Li  3p  and  3d  diffuse  functions  to 
describe  these  states,  as  mentioned  previously.  It  would  appear  then  that  even  the 
smallest  basis  set  used  provides  an  adequate  description  of  the  valence  space  and  that 
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subsequent  expansion  of  the  Li  basis  set  to  include  s,  p  and  d  Rydberg  functions  should 
yield  a  theoretical  spectrum  in  good  agreement  with  experiment.  However,  this  is  not  the 
case.  Basis  sets  2  and  3  both  contain  the  exact  same  diffuse  Gaussian  function  on  Li  to 
describe  the  Rydberg  3s  and  4s  states.  The  only  difference  between  these  two  basis  sets 
lies  in  the  degree  to  which  the  Gaussian  s-type  valence  basis  functions  are  contracted. 
Basis  set  2  contains  six  primitive  valence  Gaussian  s-type  functions  contracted  to  two, 
while  basis  set  3  contains  the  same  Gaussian  s-type  valence  basis  functions  contracted  to 
three.  This  extra  degree  of  freedom  in  the  basis  set  is  the  reason  that  only  the  3s  Rydberg 
state  appears  in  the  energy  spectrum  corresponding  to  basis  set  2  whereas  both  the  3s 
and  4s  Rydberg  states  appear  in  the  one  corresponding  to  basis  set  3.  Furthermore,  the 
energy  of  the  Li  3s  Rydberg  state  (state  14)  corresponding  to  basis  set  2  is  0.5  eV  higher 
in  energy  than  the  experimental  value  of  3.37  eVj  6  The  energy  of  the  Li  3s  Rydberg 
state  corresponding  to  basis  set  3  (also  state  14),  on  the  other  hand,  is  only  0.04  eV  lower 
than  the  experimental  value.  The  fact  that  the  energy  of  this  state  is  closer  by  a  factor  of 
ten  to  the  experimental  value  using  basis  set  3  is  a  surprising  result  given  that  both  basis 
sets  contain  the  same  Gaussian  function  to  describe  the  Li  3s  and  4s  Rydberg  states.  It  is 
clear  then  that  basis  set  2  is  too  small  to  provide  an  accurate  representation  of  the  Li  Be 
energy  spectrum,  despite  the  apparent  accuracy  of  the  resulting  Li  2P  valence  state 
energies.  Because  of  this  inadequacy,  the  diffuse  s-type  Gaussian  function  used  to 
represent  the  Li  s  Rydberg  states  is  not  only  describing  these  states,  but  it  is  also 

contributing  to  the  valence  space  in  the  case  of  basis  set  2.  The  error  in  the  energy  of 

state  14  corresponding  to  basis  set  2  is  therefore  a  result  of  the  'double  duty’  performed 
by  this  function. 

Relaxation  of  the  contraction  of  the  Li  s  basis  from  double  zeta  to  triple  zeta 
not  only  results  in  the  theoretical  energy  of  the  Rydberg  state  being  ten  times  more 
accurate  relative  to  the  experimental  value,  but  it  also  results  in  the  appearance  of  the  Li 
4s  Rydberg  state,  the  energy  of  which  is  0.02  eV  lower  than  the  experimental  value  of 
4.34  eV.1^  An  adequate  description  of  the  Rydberg  states  of  the  LiBe  system  is  pivotal  to 
the  calculation  of  an  accurate  electronic  spectrum,  a  fact  exemplified  by  the  large  errors 

in  the  energies  of  these  states  resulting  from  the  use  of  a  small  basis  set.  However,  the 

accurate  representation  of  these  states  is  affected  not  only  by  the  quality  of  the  Rydberg 
functions  incorporated  into  the  basis  set,  but  also  by  that  of  the  valence  basis  functions 
themselves.  A  good  description  of  the  electronic  states  of  the  system  requires  both 
inclusion  of  accurate  Rydberg  functions  in  the  basis  set  as  well  as  augmentation  of  the 
valence  basis  functions. 

Tabel  V  shows  that  the  energies  of  the  Li  Rydberg  3s  and  4s  states 
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corresponding  to  basis  set  4  contain  only  a  0.005  and  0.01  eV  deviation,  respectively,  from 
the  experimental  values.^  The  closer  agreement  to  experiment  results  from  the  use  of  a 
a  Slater-type  orbital  formed  from  an  expansion  of  six  Gaussian  orbitals  (STO-6G)  to 
describe  these  Rydberg  states  as  opposed  to  the  single  Gaussian  orbital  contained  in  basis 
sets  2  and  3.  The  calculated  spectrum  further  improves  when  STO-6G  representations  of 
3p  and  3d  Rydberg  orbitals  are  added  to  the  basis  set  (basis  set  5). 

Although  no  mention  has  been  made  concerning  the  adequacy  of  the  basis  set 
corresponding  to  Be,  the  same  conclusions  apply.  Be  does  not  possess  Rydberg  states  in 
the  energy  region  appearing  in  Tables  III,  V  and  VI.  However,  the  double  excitation 
states  appearing  in  Table  VI  arc  a  combination  of  both  Li  and  Be  single  excitations,  the 
energies  of  which  are  consequently  affected  by  the  quality  of  the  Be  valence  basis  set. 

The  fact  that  the  double-zeta  Li  basis  set  provides  such  a  poor  description  of  the  Li  Be 
system  leads  indirectly  to  the  conclusion  that  the  Be  basis  set  also  needs  to  be  augmented, 
and  Be  single  excitation  states,  the  energies  of  which  are  influenced  by  the  quality  of  the 
Be  valence  basis  set,  are  thereby  affected.  The  effects  of  basis  set  on  the  single  excitation 
states  of  Be  can  be  seen  by  comparing  the  energies  of  these  states  (states  5-13)  in  Tables 
III  and  VI.  The  former  table  contains  energy  values  calculated  using  basis  set  6  of  Table 
I,  while  the  values  appearing  in  the  latter  were  calculated  using  basis  set  5  of  Table  I. 
Augmenting  the  Be  valence  basis  set  from  2s2pld  to  4s3p2d  results  in  a  0.01  eV  lowering. 

As  discussed  previously,  the  Rydberg  states  arc  sensitive  to  the  quality  of  the 
valence  basis  set,  and,  as  already  observed,  it  is  possible  to  calculate  an  electronic 
spectrum  having  seemingly  accurate  energy  values  for  the  valence  states  even  with  a 
basis  set  that  is  not  of  good  quality.  All  of  the  calculations  discussed  to  this  point  were 
accomplished  at  the  dissociative  limit  of  40  bohr.  Of  course,  as  the  intemuclear  distance  is 
decreased,  the  interaction  between  centers  increases.  The  basis  set  employed,  therefore, 
must  be  balanced  enough  to  provide  a  comparable  description  of  both  centers  and  large 
enough  to  allow  the  flexibility  necessary  to  yield  accurate  electronic  state  energies.  Basis 
set  5  does  not  meet  these  criteria  because  it  is  of  triple  zeta  quality  with  respect  to  the  Li  s 
valence  states  but  only  double  zeta  with  respect  to  the  Be  s  ones.  It  was  discussed 
previously  that  a  basis  set  of  double  zeta  quality  is  too  small  to  provide  an  accurate 
description  of  the  system  at  the  dissociative  limit.  It  is  therefore  clear  that  this  type  of 
basis  set  cannot  provide  an  accurate  description  of  the  LiBe  spectrum  at  shorter 
intemuclear  distances.  The  fact  that  the  calculation  appearing  in  Table  III  and  the  one 
corresponding  to  basis  set  5  in  Table  V  yield  energies  in  such  close  agreement  to  each 
other  and  to  the  atomic  experimental  values  despite  the  poor  quality  of  the  basis  set 
employed  underlines  the  importance  of  ascertaining  that  a  proper  basis  set  has  been 
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chosen  before  proceeding  with  further  calculations,  especially  at  internuclear  distances 
near  equilibrium.  The  full  spin-orbit  procedure  used  in  the  present  calculations  is  ideal 
for  this  type  of  study  since  the  only  sources  of  error  possible  in  a  full  Cl  calculation  of 
this  type  arise  from  either  the  effective  core  potentials  and/or  the  basis  set.  If  the 
accuracy  of  the  core  potentials  is  established,  then  the  main  source  of  error  lies  in  the 
basis  set. 


C.  Configuration  Selection 

All  of  the  results  reported  so  far  have  been  calculated  using  the  full  spin- 
orbit  Cl  method,  A  typical  calculation  employing  the  highest  quality  valence  basis  set 
reported  here  and  including  one  STO-6G  orbital  for  each  of  the  four  Rydberg  states  of  Li 
(3s,  3p,  3d  and  4s)  required  approximately  25  cpu  hours  of  Cray  Y-MP  time  for  one 
internuclear  distance.  Since  the  amount  of  time  required  to  complete  a  calculation  is 
roughly  proportional  to  the  square  of  the  number  of  determinants  (which  are,  in  turn, 
proportional  to  the  number  of  orbitals),  the  size  of  the  basis  set  is  pivotal  in  determining 
the  tractability  of  a  given  full  SOCI  calculation.  For  low-lying  stales  of  the  LiBe  system, 
only  the  Li  center  requires  the  inclusion  of  Rydberg  functions  in  the  basis  and  for  this 
case,  only  four  Rydberg  states  lie  within  the  energy  range  of  interest.  Extension  of  the 
SOCI  procedure  to  systems  requiring  a  larger  valence  basis  set  and/or  the  inclusion  of  a 
greater  number  of  Rydberg  functions  to  describe  these  states  for  one  or  both  centers 
may  lead  to  unmanageable  calculations  at  the  full  Cl  level.  However,  the  systems  may 
still  be  treated  at  this  level  of  theory  providing  that  a  set  of  configurations  smaller  than 
that  of  the  full  Cl  may  be  selected  without  the  introduction  of  significant  error. 

Table  VI  shows  results  of  two  configuration  selection  schemes  -  CSS1  and  CSS2 
-  used  in  the  SOCI  calculations  for  the  LiBe  system  using  the  most  accurate  basis  set 
reported  here  (basis  set  6).  For  each  case,  the  results  are  compared  to  the  full  SOCI  to 
determine  the  accuracy  of  the  configuration  selection  scheme  employed.  The  full  SOCI 
calculations  required  9102  double  group  spatial  configurations.  CSS1  contains  all 
configurations  involving  single,  double  and  triple  excitations  among  the  valence  orbitals 
and,  in  addition,  only  those  configurations  which  correspond  to  Rydberg  states.  A  total  of 
4072  configurations  were  included.  Results  corresponding  to  thirteen  slates  are  reported. 
(The  higher-lying  states  arc  double  excitations  involving  Li  and  Be  valence  orbitals.)  No 
Rydberg  states  appeared  using  this  configuration  selection  scheme.  It  is  also  interesting 
to  note  that  the  number  of  configurations  treated  was  not  flexible  enough  to  lead  to  a 
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near  triply  degenerate  set  of  Li  2P  states.  The  spin-orbit  splitting  between  the  2?\/2  and 
2P3/2  states  is  overestimated  because  of  imbalance  in  the  number  of  configurations 
which  contribute  to  these  states.  Certain  configurations  that  involve  excitations  among 
valence  and  Rydberg  states  obviously  contribute  to  the  2P3/2  states.  Since  these  were  not 
included,  the  energy  of  these  orbitals  was  calculated  to  be  higher  than  the  actual  values, 
breaking  the  degeneracy. 

CSS2  included  all  of  the  configurations  of  CSS1,  as  well  as  all  single  excitations 
with  respect  to  each  Rydberg  configuration,  resulting  in  a  total  of  7407  configurations. 
The  results  agree  closely  with  those  of  the  full  SOCI,  the  largest  error  being  0.  00004  eV. 

In  this  case,  the  Li  2Pl/2  and  2P3/2  valence  states  were  calculated  to  be  nearly 

degenerate.  It  is  therefore  apparent  that  the  Rydberg  states  of  the  Li  atom  treated  in  this 
work  arc  truly  SCF-like.  The  radial  extent  of  these  Rydberg  orbitals  is  such  that  an 
electron  occupying  them  essentially  'sees'  a  positive  ion  core,  i.e.,  occupation  of  these 

orbitals  results  in  a  hydrogenic  system.  Inclusion  of  all  configurations  involving  single 
excitations  with  respect  to  the  Rydberg  states  is  therefore  needed  to  allow  for  ion  core 
molecular  orbital  relaxation  effects  that  arise  from  the  hydrogenic  behavior,  since  any 

configuration  space  smaller  than  the  full  Cl  results  in  the  wave  function  being 
dependent  on  the  precise  forms  of  the  orbitals.  Selectively  reducing  the  configurations 
in  the  SOCI  from  9102  to  7407  resulted  in  a  savings  of  approximately  twenty  minutes  of 
Cray  Y-MP  epu  time.  However,  since  the  time  required  to  complete  a  calculation  scales  as 
the  square  of  the  number  of  determinants  generated,  the  savings  will  prove  to  be 
significant  for  a  system  containing  a  large  number  of  Rydberg  states. 


Summary  and  Conclusions 

Ab  initio  spin-orbit  full  SOCI  calculations  on  LiBc  in  the  context  of  relativistic 
effective  core  potentials  predict  that  the  system  is  weakly  bound,  with  a  dissociation 
energy  of  0.25  cV  corresponding  to  an  equilibrium  intcmuclcar  separation  of  4.95  bohr. 
A  total  of  29  states  were  calculated  -  25  for  the  neutral  system  and  4  for  LiBe+.  Of  these, 
eighteen  appear  in  Figure  I,  while  twelve  were  spin-orbit  components  of  those  already 
plotted  and  thus  arc  not  resolved  on  the  energy  scale  shown.  The  calculated  energies  of 
all  the  atomic  states  arc  in  excellent  agreement  with  experiment,  the  largest  error  being 
0.05  cV  for  both  the  neutral  2D  +  !S  and  ionic  *S  +  !S  states.  All  of  the  excited  states  of 
neutral  LiBe  exhibit  more  bonding  than  the  ground  state.  The  energy  range  is  0.30-2.03 
cV  and  their  equilibrium  intcmuclcar  separations  range  from  4-7  bohr. 
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The  calculated  atomic  spin-orbit  splittings  are  in  excellent  agreement  with 
experiment  for  both  the  neutral  and  ionic  systems,  the  largest  error  being  O.lmeV.  The 

largest  splitting  occurred  between  the  neutral  2S  +  2Pj  and  2S  +  ^P2  and  the  ionic  +  3Pl 

and  *S  +  3P2  states.  The  rest  of  the  atomic  states  exhibit  spin-orbit  splittings  in  the  range 
0.0-0. 1  mcV.  Molecular  spin-orbit  splittings  for  the  neutral  system  range  from  0.0-0.08 
meV,  with  the  largest  splitting  occurring  between  the  components  of  the  A2n  state.  The 
largest  splitting  is  0.22  meV  between  the  components  of  the  B^n  state  of  the  cation. 

The  lower-lying  Rydberg  and  upper-lying  valence  states  were  found  to 
interact  and  affect  the  shape  of  the  reported  curves.  It  was  therefore  necessary  to 
:  ".elude  Rydberg  functions  in  the  basis  set  when  studying  the  upper  valence  region  of 
the  spectrum.  The  Rydberg  states  were  found  to  be  excellently  represented  by  STO-6G 
functions.  The  overall  accurate  representation  required  both  the  inclusion  of  accurate 
Rydberg  functions  as  well  as  augmentation  of  the  valence  basis  functions.  The  most 
accurate  basis  set  was  composed  of  4s,  3p  and  2d  functions  on  each  center  plus  STO-6G 

functions  to  represent  the  3s,  3p,  3d  and  4s  Rydberg  states  of  Li.  The  largest  error 

between  the  calculated  and  experimental  atomic  energies  using  this  basis  set  was  0.05  eV 
for  the  neutral  2D  (3d)  +  lS  state  and  0.18  eV  for  the  ionic  lS  +  JP  state. 

A  configuration  selection  scheme  (CSS)  comprised  of  a  total  of  7407 
configurations  corresponding  to  the  configurations  of  the  full  valence  Cl,  all  Rydberg 
states  configurations  and  all  configurations  corresponding  to  single  excitations  with 
respect  to  each  Rydberg  slate  was  used  to  calculate  the  atomic  states  of  LiBe.  These  were 
in  excellent  agreement  with  those  of  the  full  Cl,  indicating  that  the  Rydberg  states  of  the 
Li  atom  are  truly  SCF-like.  Results  obtained  using  this  type  of  CSS  are  comparable  in 
accuracy  to  those  of  the  full  Cl  while  reducing  the  computational  complexity.  The  SOCI 
procedure  used  in  this  study  has  been  shown  to  yield  highly  accurate  results  for  systems 

where  a  full  Cl  or  use  of  a  comparable  CSS  is  feasible. 
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Table  I. 

Basis  sets  of  Gaussian-type  functions.2 

Basis  Set 

Li 

Be 

1 

(6s5pld)/[2s2pld]b 

(6s5pld)/[2s2pld]b 

2 

(7s5pld)/[3s2pld]c 

(6s5pld)/[2s2pld]b 

3 

(7s5pld)/[4s2pld]d 

(6s5pld)/[2s2pld]b 

4 

( 1 2s5p  1  d)/[4s2p  1  d  ]e 

(6s5pld)/[2s2pld]b 

5 

(18sllp7d)/[5s3p2d]f 

(6s5pld)/[2s2pld]b 

6 

(18sllp8d)/[6s4p3d]g 

(6s5p2d)/[4s3p2d]h 

a  Parentheses  (brackets)  designate  primitive  (contracted)  basis  sets. 
b  Valence  basis  set  plus  one  polarization  function. 
c  Basis  set  1  plus  one  diffuse  s-type  function. 

d  Basis  set  2  whereby  one  valence  s-type  function  has  been  uncontracted. 
c  Basis  set  3  whereby  the  diffuse  s-type  function  has  been  replaced  by  one  STO-6G  4s  Rydberg 
function. 

^  Basis  set  4  plus  three  additional  STO-6G  Rydberg  functions,  one  each  to  describe  the  3s, 

3p,  and  3d  Rydberg  atomic  orbitals, 

8  Basis  set  5  whereby  one  additional  polarization  function  has  been  added  and  the  valence 
basis  functions  have  been  uncontracted  from  [3s2p]  to  [4s3p]. 
b  Basis  set  1  whereby  one  additional  polarization  function  has  been  added  and  the  valence 
basis  functions  have  been  uncontracted  from  [2s2p]  to  [4s3p], 
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Table  II.  Exponents  and  coefficients  of  the  STO-6G  Gaussian  expansions  of 

Slater-type  orbitals  representing  the  Rydberg  states  of  LiBe. 

Orbital  Exponent  (Qa 

Exponent  (ot)b 

Coefficient0 

3s  0.14825 

0.40480 

-0.00974 

0.11656 

-0.07266 

0.04585 

-0.17162 

0.02096 

0.12898 

0.01043 

0.72886 

0.00535 

0.30133 

3  p  0.11436 

0.31226 

-0.00810 

0.08991 

-0.01715 

0.03535 

0.07370 

0.01616 

0.39651 

0.00804 

0.49781 

0.00413 

0.11748 

3d  0.11110 

0.30335 

0.00663 

0.08735 

0.05958 

0.03435 

0.24019 

0.01570 

0.46481 

0.00781 

0.34341 

0.00401 

0.05389 

4s  0.07712 

0.10529 

0.00378 

0.03388 

-0.05586 

0.01447 

-0.31930 

0.00722 

-0.02765 

0.00390 

0.90492 

0.00217 

0.34063 

a  Slater  exponent, 
b  Gaussian  exponent, 
c  Ref.  14, 


Tabic  III. 

Energies  of  LiBc  and  LiBe+  at  the 

dissociative  limit. 

i 

Root 

Atomic'3 

Full  SOCIc 

Expt.d 

!errorie 

No. 

States 

1 

2S  +  1S 

0.0000 

2 

2Pl/2  +  !S 

1.8439 

1.8478 

0.0039 

3-4 

2P3/2  +  !S 

1.8439 

1.8478 

0.0039 

5 

2S  +  3p0 

2.7666 

2.7250 

0.0416 

6-8 

2S  +  3Pi 

2.7667 

2.7251 

0.0416 

9-13 

2S  +  3P2 

2.7669 

2.7254 

0.0415 

14 

2S  (3s)  +  lS 

3.3325 

3.3730 

0.0405 

15 

2P 1/2  (3p)  +  JS 

3.7974 

3.8342 

0.0368 

16-17 

2P3/2  (3p)  +  IS 

3.7975 

3.8342 

0.0367 

18-19 

2D3/2  (3d)  +  !S 

3.8288 

3.8785 

0.0497 

20-22 

2D5/2  (3d)  +  IS 

3.8288 

3.8785 

0.0497 

23 

2S  (4s)  +  !S 

4.3153 

4.3408 

0.0454 

1 

!S  +  JS 

5.3408 

5.3879 

0.0471 

2 

*S  +  3p0 

8.1075 

8.1129 

0.0054 

3-5 

JS  +  3Pi 

8.1075 

8.1130 

0.0055 

6-10 

IS  +  3P2 

8.1077 

8.1133 

0.0056 

11 

2S  +  2S 

9.3383 

9.3224 

0.0159 

12-14 

!S  +  !Pi 

10.8476 

10.6689 

0.1787 

15 

2P  i/2  +  2S 

11.1822 

11.1702 

0.0120 

16-17 

2p3/2  +  2S 

11.1822 

11.1702 

0.0120 

a  Energies 

are  in  eV.  The  dissociative 

limits  are  40 

bohr  for  LiBe  and 

10^  bohr  for  LiBe+. 

b  First  and  second  entries  correspond 

to  Li  and  Be 

respectively. 

c  Basis  set  6  of  Ta'Je  I. 

d  Ref.  15. 

e  Absolute 

values  of  the  energy  differences  between 

entries  in  columns 

3  and  4. 
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Tabic  IV. 

Molecular  constants 

of  electronic  states 

of  LiBe 

and  LiBe+. 

State 

Rea 

Deb 

tocc 

coexcc 

Bec 

aec 

G0C 

Tec 

X  2£f 

4.95 

0.25 

298 

9.8 

0.622 

0.021 

146 

a  2n 

4.27 

1.50 

491 

3.6 

0.838 

0.0092 

245 

4768 

b  2i+ 

4.91 

0.89 

345 

2.0 

0.633 

0.0077 

172 

9745 

c2n 

5.77 

0.64 

191 

-0.25 

0.459 

0.0075 

96 

19138 

a4n 

4.74 

0.60 

360 

4.9 

0.680 

0.012 

179 

19525 

d  2r* 

6.73 

0.30 

160 

7.1 

0.337 

-0.046 

85 

21884 

b  4i- 

3.99 

2.03 

495 

3.4 

0.961 

0.012 

247 

22836 

2L+(3s) 

5.23 

0.59 

383 

9.2 

0.558 

0.011 

188 

24173 

2L+(3p) 

5.10 

0.68 

257 

2.2 

0.587 

0.016 

129 

27157 

2n  (3p) 

4.76 

0.63 

356 

6.7 

0.673 

0.013 

176 

27583 

2A  (3d) 

4.89 

0.53 

300 

0.35 

0.638 

0.0081 

151 

28654 

2n  (3d) 

5.60 

0.47 

211 

3.1 

0.487 

-0.034 

106 

29113 

2I+(3d) 

4.81 

0.38 

357 

5.0 

0.660 

0.013 

177 

29805 

2£+(4s) 

5.03 

0.52 

302 

6.0 

0.603 

0.011 

149 

32426 

4.98 

0.56 

311 

4.8 

0.616 

0.013 

154 

-  . 

a  3I+ 

5.58 

0.91 

274 

2.2 

0.490 

0.006 

137 

19476 

b  3n 

5.09 

0.19 

235 

7.0 

0.590 

0.021 

116 

25332 

a  bohr. 
b  eV. 
c  cm’*. 
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Tabic  V.  Excitation  energies  (eV)  showing  effects  of  basis  set  on  the  Rydberg  states  of  Li  Be 
at  the  dissociative  limit  (40  bohr)  using  full  SOCI. 


Root2 

Basis  Set  1 

Basis  Set  2 

Basis  Set  3 

Basis  Set  4 

Basis  Set  5 

1 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

2 

1.8389 

1.8389 

1.8417 

1.8417 

1.8419 

3 

1.8389 

1.8389 

1.8417 

1.8417 

1.8419 

4 

1.8389 

1.8389 

1.8417 

1.8417 

1.8419 

5 

2.7761 

2.7762 

2.7761 

2.7764 

2.7764 

6 

2.7765 

2.7765 

2.7765 

2.7765 

2.7765 

7 

2.7765 

2.7765 

2.7765 

2.7765 

2.7765 

8 

2.7765 

2.7765 

2.7765 

2.7765 

2.7765 

9 

2.7767 

2.7767 

2.7767 

2.7767 

2.7767 

10 

2.7767 

2.7767 

2.7767 

2.7767 

2.7767 

11 

2.7767 

2.7767 

2.7767 

2.7767 

2.7767 

12 

2.7767 

2.7767 

2.7767 

2.7767 

2.7767 

13 

2.7767 

2.7767 

2.7767 

2.7767 

2.7767 

14 

4.3889 

3.8779 

3.4179 

3.3682 

3.3354 

15 

4.3889 

4.3889 

4.3169 

4.3508 

3.7987 

16 

4.3889 

4.3889 

4.3916 

4.3916 

3.7987 

17 

— 

4.3899 

4.3916 

4.3917 

3.7987 

18 

— 

— 

4.3917 

4.3917 

3.8268 

19 

— 

-- 

— 

3.8268 

20 

— 

— 

— 

3.8268 

21 

— 

— 

— 

3.8268 

22 

— 

— 

— 

3.8268 

23 

— 

-- 

— 

4.3153 

a  Refer  to  Table  III  for  state  designations. 


Tabic  VI. 

LiBe  excitation  energies  (eV)  corresponding  to  SOCI  wavefunctions 
defined  using  two  configuration  selection  schemes  (CSS).a 

Rootb 

Atomic  State 

SOCI/CSS1c 

SOCI/CSS2d 

Full  SOCIe 

1 

2S  +  !S 

0.0000 

0.0000 

0.0000 

2 

2Pl/2  +  lS 

1.8419 

1.8419 

1.8419 

3-4 

2P3/2  +  *S 

2.0595 

1.8419 

1.8419 

5 

2S  +  3Po 

2.7765 

2.7764 

2.7765 

6-8 

2S  +  3Pi 

2.7765 

2.7765 

2.7765 

9-13 

2S  +  3P2 

2.7767 

2.7767 

2.7767 

14 

2S  (3s)  +  !S 

4.5655 

3.3354 

3.3354 

15 

2P 1/2  (3p)  +  lS 

— 

3.7987 

3.7987 

16-17 

2P3/2  (3p)  +  !S 

— 

3.7987 

3.7987 

18-19 

2Ds/2  (3d)  +  IS 

— 

3.8269 

3.8268 

20-22 

2Dj/2  (3d)  +  IS 

— 

3.8269 

3.8268 

23 

2S  (4s)  +  IS 

— 

4.3154 

4.3153 

24 

2P 1/2  +  3P0 

4.6183 

4.6183 

4.6183 

25-27 

2Pl/2  +  3Pl 

4.6183 

4.6183 

4.6183 

50 

2P3/2  +  3P2 

4.6523 

4.6187 

4.6187 

a  Roots  24  through  50  are  double  excitation  states  of  LiBe.  Only  the  first  four  and  the 
last  of  these  states  are  shown. 


b  All  calculations  were  accomplished  at  40  bohr  using  a  triple  zeta  (double  zeta)  s,  double 
zeta  (double  zeta)  p  and  single  zeta  (single  zeta)  d  valence  basis  set  on  Li  (Be)  plus  3s, 
3p,  3d  and  4s  STO-6G  functions  to  represent  Rydberg  states  of  Li  (Basis  set  5  in  Table  I). 
c  4072  configurations.  Includes  full  valence  SOCI  plus  Rydberg  state  configurations, 
d  7407  configurations.  Includes  full  valence  SOCI,  Rydberg  state  configurations  and  all 
single  excitations  relative  to  Rydberg  states  where  no  more  than  one  electron  occupies  a 
particular  Rydberg  orbital  for  any  given  configuration. 
e  9102  configurations. 
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Figure  captions 

Figure  1.  Potential  energy  curves  of  valence  and  Rydberg  states  of  LiBe  and  of  three 
states  of  LiBe+.  Asymptotic  atomic  states  are  Li  and  Be,  respectively. 

Figure  2.  Selected  potential  energy  curves  of  LiBe  plotted  on  an  expanded  energy  scale 
relative  to  that  of  Fig.  1. 
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INTRODUCTION 

Core  polarization  is  an  effect  due  to  the  correlation  between  electrons  in  the  valence 
space  and  those  in  the  core  region.  As  a  result  of  this  correlation,  valence  electrons  are 
affected  by  the  instantaneous  positions  of  the  core  electrons,  and  the  simple  concept  of  core¬ 
valence  separability  breaks  down  as  the  magnitude  of  this  correlation  phenomenon  increases. 
Core  polarization  is  greatest  for  systems  containing  a  small  number  of  valence  electrons  and 
having  cores  that  are  polarized-easily. 

Born  and  Heisenberg  presented  a  treatment  of  the  core  polarization  effect  for  atones 
containing  a  single  valence  electron.  1  Mayer  and  Mayer  ^  and  van  Vleck  and  White! aw? 
later  presented  improved  core/ valence  polarization  potentials.  However,  the  potentials  were 
not  properly  behaved  at  short  electron-nucleus  distances,  leading  Bierman4.  Bierman  and 
Lubeck^  and  Bales'^  to  introduce  empirical  cut-off  functions  to  address  the  problem.  Since 
these  developments  a  large  number  of  researchers  have  further  contributed  by  discussing  - 
among  many  related  topics  -  two-electron  interactions,  non-adiabalic  corrections,  and 
theoretical  potentials  using  non-empirical  functions,  as  well  as  presenting  a  host  of 
applications.^  Muller  et  al.7  have  presented  a  procedure  whereby  the  core  polarization  is 
described  by  a  semi-empirical  potential  used  in  all-electron  SCF  and  valence  Cl  calculations. 


Stoll  ct  al.  have  used  their  method  in  pseudopotential  calculations .8  Christiansen^  later 
paralleled  the  procedure  of  Muller  et  al.  to  develop  an  ab  initio  core/valence  polarization 
potential  (CVPP)  in  the  context  of  the  adiabatic  field  approximation  i*3*1 1  and  relativistic 
effective  potential  (REP)  theory.  '2 

The  latter  CVPP  has  several  advantages.  First,  the  procedure  is  non-empirical. 
Second,  relativistic  effects  are  incorporated.  Third,  the  valence  space  is  reduced  to  include 
the  smallest  possible  number  of  electrons.  The  last  two  points,  which  are  extremely 
important,  result  from  derivation  of  the  CVPP  within  the  REP  approximation. 

REP-based  methods  rely  on  the  principle  of  core/valence  separability.  Consequently,  a 
choice  must  be  made  regarding  the  number  of  electrons  to  be  treated  explicitly.  In  the 
procedure  developed  by  Lee  et  al. 1 this  choice  is  arbitrary.  For  example,  the  most  practical 
REP  for  cesium  would  correspond  to  a  one-electron  valence  space.  The  same  atom  may  also 
be  described  by  9-electron  and  19-electron  valence  spaces,  where  the  former  REP  contains 
the  5s  and  5p  Cs  electrons  in  the  valence  space  and  the  latter  contains  a  valence  space 
comprised  of  the  5s,  5p  and  4 d  electrons.  A  one-electron-valence-space  REP  would  allow 
the  minimum  number  of  electrons  to  be  treated  explicitly  and,  therefore,  would  render 
calculations  at.  say.  the  full  configuration  interaction  (Cl)  level  possible  for  clusters  of  two  or 
more  Cs  atoms.  However,  such  an  REP  yields  inaccurate  excitation  and  ionization  energies 
for  Cs  because  of  the  large  effects  of  core/valencc  polarization  in  this  atom.13  To  improve 
accuracy,  a  9-electron-valcnce-space  REP  may  be  used,  although  this  would  obviously 
render  full  Cl  calculations,  or  even  procedures  that  incorporate  modest  levels  of  electron 
correlation,  intractable.  There  is  a  second  disadvantage  to  using  a  larger  valence  space  REP. 
Because  the  REPs  arc  generated  using  Dirac-Fock  theory,  relativistic  effects  are  incorporated 
into  REP-based  calculations  via  the  core  electrons.  Relativistic  effects  in  the  valence  region, 
which  is  treated  non-relalivistically,  result  from  the  propagation  of  these  effects  from  the  core 
to  the  valence  electrons.1''  Thus,  a  reduction  in  the  number  of  core  electrons  results  in  a 
loss  of  relativistic  effects,  leading  to  decreased  accuracy  for  calculations  on  heavy-element 
systems  where  such  effects  are  increasingly  important.  Christiansen  showed  that  a  high 
level  of  accuracy  may  be  achieved  with  the  minimum-valence-electron  REP  if  a  CVPP  is 
included  in  the  calculations.^ 

This  paper  extends  the  method  presented  by  Christiansen  to  an  REP-based  CVPP 
derived  in  j-j  coupling.  Here,  the  relativistic  CVPP  (RCVPP)  is  calculated  using  individual 
orbital  contributions  rather  total  energies,  making  it  possible  to  see  the  effect  of  each  orbital 
on  the  total  potential. 


METHOD 

Since  the  formal  derivation  of  the  REPs  is  given  in  Refs.  12  and  14,  only  those 
relationships  relevant  to  the  RCVPP  development  are  given  here.  The  REPs,  which  are 
based  on  numerical  Dirac-Fock  (DF)  wavefunciions15,  are  expressed  as  follows: 


L-l  1+1/2  j 

IjREP  =  UuRH>(r)  +Ii=oXj=il-i/2l  Xm«-j  (  UijREP(r)  -  Uu1^  (r)]  lljmxljml  (1) 

where  ihe  projection  operators  lljmxljml  dictate  that  UREP  must  be  used  in  the  context  of  a 
basis  set  of  two-component  Dirac  spinors. 16  The  REPs  can  be  recast  into  the  form  of  an 
angular  momentum  averaged  potential  that  may  be  used  in  the  context  of  standard 
nonrelativistic  electronic  structure  procedures  based  on  atomic  L-S  coupling  and  a  spin-orbit 
operator. 14  This  AREP  is  defined  by 

L-l  I 

UAREP  =  uLAREP  (r)  ♦  I1=0  [  U,AREP  (r) .  ULAREP  (r)]  j|mxlml  (2) 


where 


REP  REP 

ULAREP=  (21+1)1  [ lUM.in  +  (1+1)  Uu*,«]  (3) 

In  Eqs.  (1)  and  (2)  UuREPand  Ulrep  refer  to  the  so  called  "residual"  EPs  where  L  J  and  L 
are  normally  taken  as  one  larger  Ilian  Ihe  largest  angular  momentum  quantum  numbers  among 
the  core  electrons. 

The  resulting  spin-orbit  operator  is 

L  REP  I 

Hso  =  S  •  E)=i  [2/(21+ 1)]  AU|  (r)  Zm=.|  II m>  <!ml  1  llm>  <lmt  (4) 

where 

REP 

H|so  =  (2/(21+!))  AU|  (r)  (5) 


and 


ARE1  REP  REP 

AUl  (r)  =  U|,M/2  -  U|_|»ir2  (6) 

Although  U,jREPt  U|AREP_  an£j  H|SO  are  derived  as  numerical  operators,  consistent  with 
their  origin  from  numerical  DF  wavefunctions,  they  are  re-expressed  as  expansions  in  terms 
of  Gaussian  functions  in  the  following  form: 

M  2 

Ok  (r)  =  r-2  Zi=,  CU  •  <rki)  e-°u  r  .  Ok  (r)  =  U,jREP  u,AREP  or  H,SO  (7) 

The  RCVPP  is  defined  in  the  context  of  the  adiabatic  field  approximation  and  the  REP 
of  Eq.  (1).  The  total  core  polarization  potential  operator  is,8 

VrcvPP=  -1/2  Icdcfc  fc 


(8) 
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where  ac  is  the  polarizability  of  the  outer  core  of  atom  c  and  fc  is  the  electric  field  at  the 
nucleus  due  to  the  valence  electrons  and  the  other  atomic  cores.  If  an  adiabatic  potential  is 
defined  as 

VRCvppflJ(r)=Epoi(r)-Eion(r)  (9) 

Where  Ej0n  and  Epoi  are  the  energy  of  the  atom  in  the  field  of  a  point  charge  placed  at  r  before 
and  after,  respectively,  orbital  (self-consistent)  relaxation  calculated  using  a  relativistic 
procedure.  The  RCVPP  can  then  be  defined  as 

VRcvpp(r)=  C  (rci)  Vrcvpp"11^)  (10) 

1"  Eq.  (10)  the  fraction  of  core  electron  density  within  the  valence  radius  is 

Nc  rd 

C(rci)=l/Ncrk=|/ l<K(r)|2dr  (11) 

0 

where  4>k  are  the  Nc  normalized  radial  numerical  DF  core  orbitals.9 
The  effective  field  at  the  core  of  atom  c  is  then 

fc  (rci)  =  Oci/rci)  [  2 C  (rd)  (r) /ac)1'2  (12) 

which  leads  to  the  following  form: 

Vrcvpp=  ^2  2,  ac  [Si  fc(rci)2  +  Sj  fc(rci)  •  fc (nj) 

-2SiSc.fc(rQ)-(Rcc-/Rc<:,3)  Z  g 

+  Sc'.c”  ( Red '  Rcc-)  f  (Rcc3  Rcc*3)  Zc  Zc"  ]  (13) 

where 

ac  =  -2r>AEpo]/q2  (14) 

and  ac  is  the  polarizability  (sec  Eq.  8),  q  the  point  charge  at  a  distance  r  from  the  nucleus  and 
AEpol  the  change  in  energy  (Vrcvpp'k])  due  to  polarization  of  the  outer  core  orbitals  by  the 
point  charge.  In  Eq.  (13)  the  terms  one  through  four  are  the  one-center  one-electron 
contributions,  one-center  two-electron  contributions,  electron-other-core  interactions,  and 
core-core  interactions,  respectively.  Zc  is  the  charge  on  core  c.  It  is  noted  that  substituting 
Eq.  12  into  Eq.  13  results  in  the  cancellation  of  ac  from  the  first  term  -  the  one-center  one- 
electron  contribution. 

Up  to  this  point  the  development  of  the  core/valence  polarization  potential  Vrcvpp  is 
the  same  as  that  in  Ref.  10  for  the  nonrelativistic  treatment.  We  will  now  show  how  the 
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Vrcvpp  can  be  derived  from  individual  orbital  contributions.  The  tot*1  *"ergies  Epot  (R)  and 
Eion  (R)  given  in  Eq.  9  can  be  expressed17  as 

Epoifr)  =  I,  Ui +  »,,)  + ZcZVr  (15) 


and 


Eion(r)  =  Ij  (E|°  +  Hjj0)  +  ZcZ7r  (16) 

where  e;  and  Hjj  are  the  orbital  energy  and  one-electron  energy,  respectively,  corresponding 
to  a  system  comprised  of  an  atom  in  the  field  of  an  electron  or  point  charge  with  Z'  =  -l.  £,° 
and  H,j°are  the  orbital  energy  and  one-electron  energy,  respectively,  corresponding  to  this 
same  system  but  before  self-consistent  orbital  relaxation  occurs;  the  orbitals  in  this  case 
correspond  to  field-free  atomic  orbitals,  i.e.  without  the  presence  of  the  point  charge.  ZcZ'/r 
is  the  interaction  potential  the  atomic  nucleus  and  the  point  charge.  Since  Z‘  corresponds  to  a 
point  charge  (electron).  ZcZ7r  becomes  -Zc/r.  Substituting  Eqs.  15  and  16  into  Eq.  9 
yields 


VRCvpp,,d(rci)i=  (Ej  +  Hjj)  -  (e;°  +  Hjj°) 


(17) 


and 


vRCVPPml(rci)  =  Zj  VRCvppad(rci)i 


(18) 


It  is  noted  that  two-electron  interactions  are  present  via  e;  and  e;0  since 

Ej  =  Hjj  +  Xj  (2Jjj  -  Kjj)  and  e,°  =  Hjj°+ Ij  (2^°  -  Kjj°)  (19) 

where  Jjj  and  Kjj  are  the  coulomb  and  exchange  integrals,  respectively.17  Further, 

Hjj0  =  Hjjalom  +  Vj,,0  (20) 

and 

Ej0  =  £iaiom+  VjjO  (2J) 

where 

Vip°  =  <  <J>,0 1  -Z^rci  1  4h°  >  (22) 


4>i0  are  atomic  orbitals  and  ej810"1  are  the  orbital  energies  corresponding  to  the  unperturbed 


atom,  i.e.,  lhe  only  difference  between  the  unrelaxed  orbital  energies  in  the  presence  of  the 
point  charge  and  the  atomic  orbital  energies  is  the  attractive  potential  between  the  negative 
point  charge  and  the  positive  nucleus. 

Substituting  Eqs.  20  and  21  into  Eq.  17  yields 
VRCvppnd(r)  =  Ii  (q  -  £iwo"'  +  Hi  -  Hjjn,0,n  -  2Vip0)  (23) 

Eq.  23,  which  defines  Vrcvpp"11  (r)  appearing  in  Eq.  10,  shows  that  VRCvpp  can  be 
calculated  from  the  orbital  energies  and  one-electron  energies  corresponding  to  relaxed 
relativistic  atomic  orbitals  in  the  field  of  a  negative  point  charge  (e;  and  Hj;,  respectively), 
and  the  orbital  energies  and  one-electron  energies  of  the  unperturbed  atom  (e;alo,n  and 
Hiimo,M. respectively)  plus  the  attractive  potential  energy  between  the  nucleus  and  the  point 
charge  (Vj,,°). 

CALCULATIONS 

The  VRCvppB(1( total)  shown  in  Figs.  1  and  2  and  the  VftCvpp(!°tal)  in  Fig.  3  were 
calculated  from  the  total  energies  Epoj  and  Ej0n  (see  Eq.  9)  using  the  method  of 
Christiansen.9  The  same  values  were  also  calculated  from  the  individual  VRCVppai1  (r)j  and 
Eqs.  18  and  23.  As  expected,  the  two  methods  yield  identical  results.  The  VRCVpp^d  (r)j 
appearing  in  Figs.  1  and  2  and  the  Vrcvpp  (r)i  appearing  in  Fig.  3  were  calculated  from 
respective  cesium  ion  energies  defined  by  Eqs.  17-22. 


Figure  1.  Total  adiabatic  potential  and  its  individual  Ij.nij  contributions  for  Cs. 
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The  Cs+  orbital  space  was  partitioned  as  follows.  The  ls-4s;  the  m=l/2  and  3/2  components 
of  the  2p-4p;  and  the  m=l/2,  3/2  and  5/2  components  of  the  3d  orbitals  defined  the  core 
space  and  were  described  by  an  REP  derived11  from  the  orbitals  of  a  singly-charged  Cs 
positive  ion  defined  as  having  18  valence  electrons.  The  4d3/2,  4ds/2.  5s.  5pi/2  and  5p3/2 
orbitals  were  treated  as  outer  core  (or  inner  valence)  and  the  6s  electron,  comprised  the 
valence  space.  A  basis  set  of  three  5s.  two  6s,  three  5p  1/2  and  5p3/2.  two  6pj/2and  two 
6p3/2.  two  4d5/2  and  two  4d3/2,  two  5d«/2  and  two  5d3/2,  and  two  4f7/2  and  two  4f5/2  Slater- 
type  functions  was  used  to  describe  these  orbitals.  The  6s  electron  was  ionized  and  the 
effects  of  polarization  by  a  negative  point  charge  on  the  4d3/2, 4d5/2.  5s  and  5pj/2  outer  core 
orbitals  of  Cs+  were  determined  by  calculation  of  Ep<,|  and  Ejon  of  Eq.  9  from  the  individual 
orbital  energies  (refer  to  Eqs.  17-22),  resulting  in  (he  generation  of  a  RCVPP  through  which 
these  effects  are  described.  Once  the  RCVPP  is  determined,  it  can  be  used  in  calculations 
where  the  4d3/2.  4dsn.  5s,  5p  1/2  and  5p3/2  electrons  of  Cs  are  considered  part  of  the  core 
and  described  by  a  one-valence-electron  (6s)  REP  to  which  the  RCVPP  is  added,  thereby 
reducing  to  one  the  number  of  valence  electrons  treated  explicitly.  This  procedure  considers 
core/valcnce  polarization  to  be  an  atomic  property.  Hence,  once  the  RCVPP  is  generated,  it 
can  be  used  to  describe  the  structure  of  Cs  in  any  bonding  situation  and  any  electronic  stales 
involving  excitation  of  the  Cs  6s  valence  electron. 

A  relativistic  diatomic  molecule  program  system  was  used  to  obtain  the  energy 
corresponding  to  each  valence  orbital.16  It  is  based  in  w-ru  coupling18  and  yields 
wavcfunclions  that  are  eigenfunctions  of  the  total  angular  momentum  squared  operator, 
making  it  possible  to  calculate  a  Vrcvpp  for  each  j  and  nij  component.  The  one-electron  and 
attractive  potential  energies  were  calculated  using  codes  developed  specifically  for  this 
purpose.  Eq.  17  was  then  used  to  obtain  each  VRcvp|>nd  (r)i-  The  set  of  C(ro)j  given  by  Eq. 
11  were  calculated  using  the  normalized  radial  numerical  D  ic-Fock  core  orbitals  $k-15 
These  same  orbitals  were  used  to  derive  the  18- valence  electron  REP  employed  here.  The 
polarizability,  ac  of  Eq.  14,  which  is  a  constant  for  sufficiently  large  r,  was  calculated  by 
taking  r  as  20  bohrs.  The  final  VRCvt>p was  calculated  using  Eqs.  23  and  10-14. 


DISCUSSION 

There  are  advantages  to  deriving  VRCvpp  from  the  relativistic  atomic  orbitals  rather  than 
from  the  total  energy.  First,  the  former  method  allows  the  contributions  of  the  VRCVpp(r)i's 
to  the  total  vncvpp10  be  studied  individually,  while  the  latter  one  does  not.  Second,  the 
calculations  are  carried  out  in  thej-j  coupling  scheme,  thereby  making  it  possible  to  determine 
contributions  of  VRcvpp  to  the  spin-orbit  splitting  energies.  Figure  1  shows  the  behavior  of 
each  Ijm  component  of  VRcvpf/“1(r).  The  five  VRcvpp"J(r)ijm shown  in  Fig.  2  were  obtained 
by  summing  over  the  m  components  of  each  Vj(CVpj/’,,(r)ij  (see  Fig.  1).  The  total  sum  of  the 
five  potentials  in  Fig.  2  and  the  potential  generated  from  the  total  energy  of  the  system 
(VRCVPPatl(<otal)],  which  are  identical,  are  also  shown  in  the  figure.  Figure  3  shows  the 
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individual  VRCvpp(r)ij  calculated  from  Eq  (10)  and  the  VRcvppntl(r)jjin  of  Fig.  2.  The 
Vrcvpp( total),  which  was  generated  from  the  individual  orbital  energies,  is  also  included. 

Since  the  complete  Vpcvppad  is  essentially  the  energy  difference  between  the  relaxed 
and  unrelaxed  orbitals  in  the  presence  of  the  point  charge  [multiplied  by  the  cut-off  function, 
of  course,  in  the  case  of  VRcvpp(r)j].  a  positive  Vrcvpp”41^);  for  a  given  orbital  indicates  that 
the  orbital  is  destabilized  by  the  charge,  while  a  negative  Vrcvpp^V);  indicates  a  stabilization 
of  the  orbital  due  to  the  presence  of  the  point  charge.  Comparing  the  VRcvpp(r)ijm  and 
V Rcvppad  (total)  in  Fig.  1,  it  is  clear  that  many  of  the  details  relating  to  orbital  behavior  in  the 
presence  of  the  point  charge  are  lost  if  Vrcvpp  is  calculated  from  the  total  energies  rather  than 
from  the  individual  orbital  energies.  This  Figure  shows  that  the  P3/23/2  and  P3/2.1/2  orbitals, 
which  arc  degenerate  at  r=“  (no  perturbation  due  to  the  point  charge),  polarize  differently. 
The  P3/2.1/2  orbital  is  more  sensitive  to  the  point  charge  than  the  P3/23/2-  The  positive  values 
of  VRcvppnt!(r)p3/2.3/2  °nd  VRcvppnd(r)p3/2. 1/2  indicate  that  the  corresponding  orbitals  are 
destabilized  by  the  point  charge.  On  the  other  hand,  the  orbitals  corresponding  to 
VRcvppad(r)dM.M.  VRCvpp^Od 5/2312  and  VRcvpp^Ods^.itt  are  all  stabilized  strongly  by 
the  field  due  to  the  point  charge  (relative  to  the  other  V|tCvppatl(r)ijm  of  Fig.  1)  since  these  all 
have  negative  values.  However,  the  differences  among  these  orbitals  are  not  large  for  r 
>1.6  bohr.  At  smaller  r.  the  behavior  of  Vrcvpp”41  (r)ijm  for  the  d  orbitals  is  mostly  dictated 
by  the  65/2.1/2  component,  which  is  not  affected  by  the  point  charge  at  r  ■<  0.8  bohr.  It  is 
noted  that  multiplication  of  each  VRcvpp!Kl(r)i  by  the  cut-orf  function  will  alter  the  behavior  of 
these  potentials  at  small  r;  this  is  discussed  further  on.  The  d3/2.i/2  orbital  is  also  stabilized 
by  polarization  and  reacts  similarly  to  the  63/23/2  orbital,  causing  the  two  corresponding 
vRCVppn<I(r)i  to  be  similar  in  shape,  although  the  63/23/2  component  is  slightly  more 
stabilized.  The  nondcgcncrale  s  and  pi/2  orbitals  also  react  similarly  to  the  point  charge,  but 
whereas  the  effect  on  the  former  is  almost  negligible,  the  latter  is  somewhat  destabilized. 

The  sum  of  the  m  components  of  each  ^RCvr^kOljm  appear  in  Fig  2.  Since  the  s  and 
pi/2  orbitals  are  degenerate,  there  is  no  difference  between  their  corresponding  Vrcvpp”11  (r)ij 
and  VRcvppad  (r)|jm.  The  65/2  and  P3/2  orbitals,  on  the  other  hand,  react  strongly  to  the 
presence  of  the  point  charge,  with  the  (1&2  orbital  being  more  affected  by  the  polarization 
process.  A  comparison  of  Figs.  1  and  2  indicates  that  the  value  of  Vrcvpp”41  for  P3/2  largely 
arises  from  the  behavior  of  the  P3/2.1/2  component.  The  maximum  at  small  r  in  the  curve 
corresponding  to  the  65/2  orbital  is  due  to  Vrcvpp041  (r)d5/23/2'  From  Fig.  2  it  is  apparent  that 
Vad(tota!),  which  is  the  sum  of  each  is  VRcvppml(r)ij,  is  closer  in  shape  to  the  VRcvpp”ll(r)d3/2 
and  VRcvppn<1(r)d5/2  potentials.  The  total  potential  indicates  that  the  overall  effect  of 
polarization  by  the  point  charge  is  to  stabilize  the  Cs  inner  valence  (outer  core)  orbitals,  while 
the  breakdown  of  VnJ(iotal)  clearly  shows  that  only  the  d  components  become  more  stable. 
Both  the  p  orbitals,  on  the  other  hand,  are  destabilized,  with  the  bulk  of  the  destabilization 
arising  from  the  polarization  of  the  P3/2  orbital. 

The  individual  VRcvpi>"4l(r)ijof  Ftg.  2 were  multiplied  by  the  cut-off  functions,  Et|s. 
10-11,  to  ensure  proper  behavior  as  r  approaches  zero  and  the  results  were  then  plotted  in 
Fig.  3,  along  with  the  total  Vrcvpp  obtained  according  to  Eqs.  18,  10  and  23.  The  resulting 


VRCVPp(lo!al)  differs  considerably  from  Vn()  (total )  of  Fig.  2.  The  shape  of  Vrcvpp  (total)  is 
now  closer  to  that  of  the  s  potential  than  to  VRcvpp(r)d3/2  and  VrcvppMjm.  as  was  the  case 
in  Fig.  2.  Multiplication  by  the  cut-off  function  also  results  in  a  smoother  ds/2  potential, 
relative  to  that  appearing  in  Fig.  2.  such  that  the  d5/2  and  d3/2 components  behave  similarly. 
Again,  the  VRcvpp(r)p3a  indicates  a  destabilization  of  its  corresponding  orbital  with  respect 
to  polarization,  although  the  instability  is  greatly  reduced  for  VRcvpp(r)p3/2  vs. 
Vncvppa<1(r)p3/2-  The  pi/2  orbital  is  now  only  negligibly  affected  by  the  point  charge, 
whereas  the  s  orbital  is  slightly  more  stable  relative  to  the  unperturbed  atom.  Again,  it  is 
clear  that  the  behavior  of  VRCVpp(totaI)  mirrors  the  behavior  of  the  VRcvpp(r)s  and.  to  a  lesser 
degree,  that  of  the  VRcvpp(r)pi/2.  However,  it  neither  indicates  the  degree  to  which  the  pj/2 
orbital  destabilizes  nor  the  d3/2  and  dsn  orbitals  stabilize  as  a  result  of  polarization.  Only  the 
Vticvpp(r)|j  were  plotted  (as  opposed  to  Vftcvpt-(*)ljm)  because  there  is  a  one-to-one 
anology  between  these  values  and  the  corresponding  Ij  terms  in  the  REP  expansion.  A 
method  for  incorporating  VRcvppfr)ijm  into  REP-based  molecular  calculations  is  currently 
under  investigation,  although  this  is  a  more  complex  problem. 

Having  developed  the  procedure  whereby  a  relativistic  Vrcvpp  is  calculated  Tor  each  Ijm 
component,  the  potentials  will  be  used  in  calculations  such  that  each  angular  momentum 
component  of  the  total  REP  has  associated  with  it  a  VrcvpKOij  °r  vRCvpKr)lji».  The  current 
procedure  involves  incorporation  of  the  total  Vrcvpp'3 and  not  its  individual  components 
(Eq.  17).  Because  the  individual  Vrcvpp  (r)j  vary  widely  from  the  total  potential,  it  is 
anticipated  that  incorporation  of  the  Ij  or  Ijm  components  of  Vrcvpp  will  improve  calculalional 
results  on  systems  where  corc/valence  polarization  effects  are  important. 
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